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1 . INTRODUCTION 


In  the  first  report  in  this  sequence  we  presented  results  of 
numerical  calculations  in  which  we  simulated  shock  propagation  in  a 
three-dimensional,  face-centered-cubic  (FCC)  lattice.  The  details  of 
the  calculation  were  not  discussed  in  order  that  the  reader  not  be  lost 
in  a welter  of  detail,  superfluous  to  understanding  the  significance  of 
the  results.  In  this  report  we  describe  more  fully  the  calculations 
undertaken  and  the  computer  program  used.  The  material  herein  is  in- 
tended primarily  for  those  who  may  wish  to  perform  similar  calculations. 

We  begin  in  Section  2 by  describing  the  lattice  used  in  the  cal- 
culations. Included  in  the  description  are  the  labeling  convention  used 
to  identify  a particular  atom,  the  method  for  finding  its  nearest 
neighbors,  the  method  for  calculating  its  potential  interaction  with 
other  atoms  in  the  lattice,  and  the  method  for  calculating  the  lattice 
constant.  In  Section  3,  we  discuss  how  initial  conditions  of  the  lattice 
(prior  to  shock  compression)  are  determined  and  how  the  shock  wave  affects 
the  assumed  periodicity  of  the  lattice.  The  differential  equations  of 
motion  satisfied  by  the  atoms  are  also  written  down.  In  Section  4,  we 
describe  the  thermodynamic  variables  (pressure,  temperature,  flow  veloc- 
ity, and  density)  calculated  in  the  program  and  explain  how  th'  r calcu- 
lation is  carried  out.  The  Appendix  discusses  the  computer  program  used 
and  provides  a listing. 


2.  DESCRIPTION  OF  THE  LATTICE 

The  model  for  the  simulations  consists  of  an  infinite,  three-dimen- 
sional, FCC  lattice  constructed  of  cells  such  as  that  shown  in  Figure  1. 
The  cube  edge,  aQ,  is  hereafter  referred  to  as  the  lattice  constant. 

Prior  to  compression,  the  crystal  is  assumed  to  obey  periodic  boundary 

conditions  in  all  three  Cartesian  directions.  Therefore,  if  a particle 

is  located  at  point  (x,y,z),  there  is  a corresponding  particle  at  the 

position  (x  t IL  a , y + n Li  , z+nLa)  where  L , L and  L are 
x o y o zo  xy  z 

positive  integers  which  define  the  periodicity  of  the  lattice  in  the 

three  directions,  and  i,  m,  n are  arbitrary  integers.  Furthermore, 

for  any  function  F which  depends  on  the  particles'  velocities  and  their 

displacements  from  equilibrium,  we  have 

F(x+  *LxaQ,  y+mLyao,  z+nLzaQ)  =•  F(x,y,z)  . (2.1) 

T.  J.D.  Powell  and  J.H.  Batteh,  "Shock  Propagation  in  the  Three- 
Dimensional  Lattice.  I.  Model  and  Results",  BRL  Report  (to 
be  published  concurrently  with  this  report). 


It  is  obvious  from  Kq.  (2.1)  that,  if  the  lattice  is  periodic,  it  is 
completely  specified  by  the  displacements  and  velocities  of  all  atoms 
within  a given  period.  We  will  henceforth  refer  to  the  period  made  up 
of  atoms  whose  equilibrium  positions  satisfy  the  relations 


0 < x < La 
x o 

0 < y < La 
1 y o 

0 < z < La 

Z O 


(2.2) 


as  the  "primary"  lattice. 

2.1  Location  of  Hquilibrium  Positions  and  Labeling  Convention. 


The  crystal  axes  will  be  assumed  to  be  coincident  with  the  coordi- 
nate system.  If  the  lattice  is  at  equilibrium,  then,  it  is  evident 
that  a plane  of  atoms  located  at  z*naQ,  where  n is  an  integer,  is  ori- 
ented with  respect  to  the  x and  y axes  as  shown  in  Figure  2a;  a plane 
located  at  z ■ (n-*-*i)ao,  on  the  other  hand,  is  oriented  as  shown  in 

Figure  2b.  For  purposes  of  calculation  it  is  most  convenient  to  label 
each  atom  by  a single  set  of  indices  (i,j,k)  and  to  do  so  we  adopt  the 
following  convention;  Atoms  in  the  first  plane,  located  at  z«0,  whose 
x and  y components  are  integral  multiples  of  a^  will  be  labeled  by  k*l; 

those  whose  x and  y components  are  half-integral  multiples  of  a will 

be  labeled  by  k=2;  atoms  in  the  second  plane,  located  at  z-a^/2,  whose  x 

components  are  half-integral  multiples  of  a and  whose  y components  are 

integral  multiples  of  afl,  will  be  labeled  by  k»3;  those  whose  x components 

are  integral  multiples  of  a and  whose  y components  are  half-integral 

multiples  of  <ZQ  are  labeled  by  k=»4.  The  process  is  then  repeated  so 

that  atoms  in  the  third  plane,  located  at  :=a^,  whose  x and  v components 

are  integral  multiples  of  are  labeled  by  k=5  and  so  forth.  Thus,  in 

general,  the  plane  in  which  particle  (i,j,k)  lies  is  given  by 


M 


(2.3) 


where  the  brackets  denote  integer  division  and  the  number  of  atoms  in 


such  a plane  of  the  primary  lattice  is  2LxLy.  The  value  of  i for  a 

particle  located  at  x»0  or  x*ao/2  is  unity  and  values  increase  along 

the  positive  x axis.  Similar  considerations  apply  for  the  index  j. 
This  convention  is  not  the  only  conceivable  one,  but  does  appear  to  be 
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I <1 

I 


Figure  2a 


Labeling  scheme  for  the  plane  of  atoms  located  at  z»na 
(For  convenience,  L -4,  L «3.)  ' 


i*4 


the  most  convenient  for  the  problem  at  harnl.  Figures  2a  and  2b  illustrate 
the  labeling  scheme  for  the  planes  located  at  z«nao  and  z-(n*S)<*0. 

respectively. 


If  we  employ  the  above  convention,  it  can  be  shown  that  the  equili- 
brium position  of  particle  (i,j,k),  normalized  by  a(l  is  given  by 

Xijk  =i_1  * S ^M0P^M0^Ck,4),4:>  j 

Yijk  aj_1  * ^ |M00(k-l,2)|  (2.4) 

ztjk -'>[¥]  vk> 

In  this  equation  MOD(I,J)  is  the  remainder  of  I divided  by  J and  the 
functions  9+  are  given  by 

0>  (k)  =1  for  k > 0 


and 


= 0 

e_  (k)  = i 


for  k < 0 
for  k < 0 


(2.5) 


= 0 for  k > 0 

Equation  (2.4)  can  also  be  inverted  to  yield  i,j,  and  k in  terms  of  the 
particle  coordinates.  The  algebra  is  rather  tedious  but  ultimately 
yields  the  result 


“ijk*1*  |W>D(2X°jk-l,2)| 
2 - v 

. = lMOD(2Y°jk-l,2)! 

k = 42°jk+l*  |MOD(2Y?jk,2)|. 


(2.6) 


2.2  Relative  Coordinates  of  Neighbors. 


For  all  reasonable  interatomic-force  models,  the  force  exerted  on  a 
particular  atom  by  the  remaining  atoms  in  the  lattice  decreases  rapidly 
with  increasing  separation  distance.  Consequently,  when  computing 
forces,  it  is  necessary  to  consider  only  atoms  in  the  immediate  vicinity 
of  the  atom  in  question.  Therefore,  given  an  atom  (i,j,k),  it  is 
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necessary  to  determine  the  values  (i' , j' , k' ) for  all  its  significant 
neighbors  in  order  to  know  which  particles  exert  an  appreciable  force 
on  the  atom. 

The  simplest  method  for  determining  an  atom'^  neighbors  is  to  note 
that,  in  the  equilibrium  lattice,  the  relative  coordinates  from  atom 
(i.j  ,k)  to  one  of  its  neighbors  are  independent  of  i,j,  and  k since  the 
lattice  possesses  translational  symmetry.  For  convenience,  then,  we 
consider  atom  (0,0,11  and  note  from  Eq.  (J.4)  that  the  relative  coor- 
dinates to  atom  (ix  ,}'  , k'  1 are  given  by 


AX°  =>  i 


/ + , ptOlH 4 + M0P Ik'  ,41,41  j 


AY°  = j'  ♦ '>  | MOD ( k7  -1,41 


(2.7) 


[t^‘]  . s [t^] 


9 (k'  ) . 


By  letting  i' ,j/ ,k7  vary  over  several  values  near  i.j,  and  k (say,  from 
-10  to  101  one  can  then  use  Eq,  (-.7)  to  determine  the  relative  coor- 
dinates to  various  neighbors. 

A short  computer  program  has  been  written  which  performs  the  above 
calculation  and  the  results  for  the  first  five  sets  of  nearest  neighbors 
are  shown  in  Table  I.  The  first  column  in  the  table  indicates  which  set 
of  neighbors  is  being  considered  and  the  last  column  the  number  of  atoms 
contained  within  the  set  (obtained  by  taking  all  possible  combinations 
of  sign  for  the  relative  coordinates).  Column  5 indicates  the  distance 
between  the  two  atoms  given  by 


1ST  = [(AX0r  ♦ (AY°) " ♦ (AZ°)  “ ] 'J 


(2.8) 


Once  i.i  and  k are  given  and  the  relative  coordinates  determined,  values 
of  i'  ,y  , and  kr  can  then  be  found  from  Eq.  (2 . t>) . 

2 . 3 Interatomic  Potential  and  Lattice  Constant. 

In  all  our  calculations  we  will  assume  that  the  atoms  in  the  lattice 
interact  through  a Morse-tvpe  potential.  Therefore,  for  a pair  of 
isolated  atoms  the  potential  energy,  normalized  by  the  dissociation 
energy,  0,  of  the  pair,  is  given  by 


-R(r/r  -1)  ] , 

e T 


(2.9) 
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In  this  expression  r is  the  separation  distance  of  the  pair,  r the 

o 

separation  at  the  minimum  of  the  potential,  and  R is  a dimensionless 
parameter  which  is  indicative  of  the  nonlinearity  of  the  potential. 


TABLE  l.  Relative  Coordinates  of  Neighbors  of  a Particular  Atom  in  the 
FCC  Lattice. 


Neighbors 

AX° 

AY° 

A 2° 

DIST 

Number 

1st  nearest 

±h 

±<S 

0 

l/S7 

12 

±H 

0 

±4 

0 

±>4 

2nd  nearest 

il 

0 

0 

l 

6 

0 

±1 

0 

0 

0 

±1 

3rd  nearest 

±1 

±s 

±»s 

Ss/2 

24 

±1 

±H 

±>s 

±1 

4th  nearest 

0 

±1 

±1 

/T~ 

12 

±1 

0 

±1 

±1 

±1 

0 

5th  nearest 

±3/2 

-*4 

0 

/s 72 

24 

±3/2 

0 

0 

±3/2 

±4 

0 

±*S 

±3/2 

±H 

±3/2 

0 

±h 

0 

±3/2 

When  the  atom  pair  is  contained  within  a lattice  the  situation  is 
more  complicated  because  the  atoms  interact  not  only  with  one  another, 
but  also  with  the  remaining  atoms  of  the  lattice.  As  noted  previously, 
for  purposes  of  computation  it  is  most  convenient  to  neglect  the 
potential  interaction  for  atoms  sufficiently  separated  that  their 
forces  of  interaction  are  negligible.  Therefore,  for  a lattice  in 
which  all  atoms  are  at  rest  in  their  equilibrium  positions,  we  have 
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for  the  potential  energy  of  atom  (i,j,k) 


n r t 

1,  / "R(|  rijk'r« 

►i  i k = 12  X e 

,3,  l.rn.n 


/rQ-l) 


2 

1 


(2.10) 


In  Eq.  (2.10)  the  prime  indicates  that  the  sum  runs  only  over  atoms 
(£,m,n)  which  are  close  enough  to  atom  (i,j,k)  to  contribute  signifi- 
cantly to  the  force  and  r?^  denotes  the  equilibrium  position  of  atom 

(i,j,k).  The  factor  one-half  is  included  so  as  not  to  count  interac- 
tions twice  when  computing  the  total  potential  of  the  lattice. 


If  we  define  a nondimensional  position  vector  according  to 

3°  - ?° 

ijk  ijk/„ 


(2.11) 


the  expression  becomes 


1/2  l 

i.  ,m,n 


-R<\ 


where 


A = CL  / T* 
o o' Io 


(2.12) 


(2.13) 


is  the  lattice  constant  normalized  by  the  equilibrium  separation  for 
two  isolated  particles.  The  data  in  Table  I can  then  be  used  to  cal- 
culate the  potential  of  particle  (i,j,k)  assuming  up  to  fifth-nearest- 
neighbor  interactions. 

In  order  to  generalize  Eq.  (2.12)  to  include  the  case  in  which 
the  atoms  are  not  in  their  equilibrium  positions,  we  define  a critical 
rad’us  Rc  such  that  interactions  between  atoms  separated  by  a distance 

greater  than  are  negligible.  The  appropriate  expression  for  the 

potential  then  becomes 


' 1/2 Z 


I ‘ 

,m,n  _ 


(2.14) 


^ijkAnJ*  Rc 


where  the  sum  extends  over  all  neighbors  which  lie  within  a sphere  of 
radius  Rc  centered  at  atom  (i,j,k).  The  zero  superscripts  have  been 

removed  from  R. and  it.  to  denote  that  these  vectors  do  not  neces- 
i j k tmn 

sarily  refer  to  the  equilibrium  positions  of  the  atoms. 
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Equation  (2. 12)  can  bo  used  to  calculate  the  nondimens ional  lattice 
constant  A in  terms  of  R by  noting  that  in  tho  equilibrium  configura- 
tion the  lattice  will  readjust  itself  so  as  to  minimize  the  potential 
with  respect  to  A^ . Therefore,  wo  must  havo 


d*.  .. 


- « l I 

l , m , n 


ljk  Cmn 


-2R(A 

e 

- 


,l*fn 


e i 

tran1 


-1) 


0 


(2.15) 


A computer  program  has  been  written  which  solves  this  equation  numeric- 
ally for  A as  a function  of  R and  tho  number  of  neighbors  taken  in 

the  sum.  Results  are  shown  in  fable  II.  The  values  of  N listed  across 
the  top  row  indicate  that  the  sum  in  F.q.  l2. 15)  included  all  atoms 
through  Nth  nearest  neighbors  and  values  of  K are  listed  in  the  first 
column. 

TABLE  II.  Values  of  the  Lattice  Constant  A as  a Function  of  R and 
Number  of  Neighbors. 


n; 

1 

' ) 

3 

4 

5 

b 

y 

8 

9 

10 

.5 

1.4142 

1.2350 

. 94lb 

.8715 

.7734 

.7472 

.6465 

.6375 

. 5905 

. 5662 

1.0 

1.4142 

1.2b29 

.9775 

.9065 

. 8070 

.7810 

.6764 

.6671 

.6186 

. 5933 

2.0 

1.4142 

1.3138 

1.0827 

1.0158 

.9204 

. 8960 

. 7875 

.7773 

. 7255 

.bR'.S 

3.0 

1.4142 

1.3522 

1.2170 

1 . 17b8 

1.1227 

1.1099 

1.0523 

l .0467 

1.0207 

1 . 0066 

4.0 

1.4142 

1.3774 

1.3172 

1.3034 

1.2898 

1.2874 

1 . 2796 

1.2791 

1.2772 

1.2764 

5.0 

1.4142 

1 . 3927 

1 . 368 1 

1 . 3640 

1.3610 

1 . 3b 06 

1 . 3596 

1 . 3596 

1 . 3594 

l . 3594 

b.O 

1.4142 

1.401b 

1.391b 

1 . 3904 

l . 3897 

l . 389b 

1 . 3895 

1 . 3895 

1 . 3895 

l . 3894 

1 | 

1.4142 

1 . 40b8 

1.4027 

1.4023 

1.4021 

1.4021 

1.4021 

1.4021 

1.4021 

1.4021 

From  the  table  it  is  evident  that  for  small  values  of  R the  po- 
tential is  extremely  long-range.  For  example,  for  R less  than  about 
3.0,  A^  has  not  converged  even  though  all  atoms  through  tenth- 
nearest  neighbors  were  employed  in  the  calculation.  Clearly  the  range 
of  the  potential  extends  beyond  this  point.  As  R increases,  however, 
the  range  of  the  potential  decreases  significantly.  For  example,  when 
R = b.  0,  Aq  has  nearly  converged  to  its  asymptotic  value  when  third- 
nearest  neighbors  are  employed  in  the  calculation.  Obviously,  the 
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higher-order  neighbors  contribute  negligibly  to  the  force.  As  R -»«, 
approaches  vT  because,  in  that  case,  the  distance  between  nearest 

neighbors  is  r and  from  Table  I 

O 


ro/a 

o 


1/  vT 


(2.16) 


or  A * n/~T. 


3.  GENERATION  OF  SHOCK  WAVE 
3 . 1 Initial  Conditions  and  Method  of  Compression. 

Before  the  lattice  described  in  the  previous  section  can  be  sub- 
jected to  shock  compression,  the  initial  displacements  and  velocities 
of  the  atoms  in  the  primary  lattice  must  be  specified.  In  this  sub- 
section, we  will  discuss  the  method  for  determining  these  initial 
conditions. 


The  most  reasonable  assumption  for  the  initial  state  of  the  lattice 
is  that  it  is  in  thermal  equilibrium  at  temperature  T.  Therefore  the 
velocities  of  the  atoms  within  the  lattice  must  be  distributed  accord- 
ing to  a Maxwellian  distribution. 


Let  us  denote  by  V. the  velocity  of  particle  (i,j,k),  normalised 

to  c/t'/m,  where  D is  the  dissociation  energy  for  an  isolated  pair  of 
atoms  and  m is  the  mass  of  the  particle.  Then  the  probability  that 

the  a**'  Cartesian  component  of  lies  between  V and  V*dV  is  given 

bv  the  distribution  function 


( fe  f 


-YV2/2 


dV. 


(3.1) 


In  Equation  (3.1),  y is  the  ratio  of  the  dissociation  energy  to  the 
thermal  energy,  viz., 


I) 

Y s kTT 


(3.2) 


where  T is  the  absolute  temperature  and  is  Boltzmann's  constant. 

Furthermore,  the  total  energy  within  the  lattice  is  a constant  of  the 
motion  and,  if  the  lattice  were  harmonic,  would  be  given  by  3nkRT 

where  n is  the  total  number  of  atoms  in  the  lattice. 


dittoes  Use  of  thi  • 

H '*<='•  so  th>t  V''»C„ y df 

- I.L  V • 0 


w.  . *.L  v • » "VUMiti 

“"ha'^,/*  '^a  notJd"Stu7s  that  „„  f; 

*•  *««  -^SMTL**«  STutT  **  ^®Parted  t 

■’««  cd ;; ;; to  e.rr»ts/;f  ^ 

^ a co™0"  Pos0;tith  t,,is  ene"‘"  - 

co  °e  consic*.  °ns  *nd  scaled 

ste"'  Kith  £.*'  «'»■ 

r V2  on  tota' 


(3.3) 


. I . 

o^callyalired  the  Part  ‘ ^ ^ ^ ’ r 

££,‘2*.  .o  oscii, 

L*'Lr'L,'<  an/„*f“«  3 ,C«  r indeed 'j1'"  distr‘buuZ^!y ' p°n 

;*•  *■*«  of .,  ,,,otta'i  °"  wcea,rer“iu,t  s*.11,*-  «2??t 

interval  aiv_n  . toms  whose  Vei  • CaI  ajcis  is  h Ce  with 

shs  *■' t2r.««rvr-*  •'o  csr-  °f 

, a„ . “"wai  --  -z:ia-6^C:;h; 

i/'te  - of  inlti. 


V2 

ijk 


„ When  a„  —■ 'Z;Z'£r.ct'T  „ 

the  lattice  j pProPriate  Set.  ^ be  Constant  j, 

assume  that  n Suhjected  * initial 

SSvSS-Sp^S^r^. 

^«S5S^2£j|:^fc- 

'»«  “'Cjan«  *«  *.o.  *.  c atl°'’ is  p"f»s«cjtr *» 


“‘'-‘ore,  tlw  ca.-  l disturi  the  ” “ttice.  ,V“‘a'i«ns 

'-re  r'f8^^  at  =.„.  „ co  “,ati°"  ia  in 

to  those  o?  ?kSerl®»  o:  se  nstruct  a sen,-  foI1°Ws: 

fa'°*s  no™.,  t ' PrlM^  '^'ce'.5  t\°  '*•««  at  t 

lip'in  Pha"C  ':’°sub/ect*d'itS'  as  ahow  l„S*<!"ent  “"**”» Wl*^ 

P Positive  2 Hi  to  s toady  rr/0  Fl*Ure  4 ..  - Pianes  of 

velocity  is  a dlrection  to^  mPression  hv  At  time  t=o 

- - 


in 





We  now  consider  the  j plane  assumed  to  be  located  in  the  first 
segment.  Before  the  shock  front  reaches  this  plane,  the  motion  of  a 
particular  atom  contained  within  it  is  the  same  as  the  motion  of  the 

corresponding  atom  in  the  j plane.  The  identical  motion  results 

because  these  atoms  had  the  same  initial  conditions  and  were  acted 

upon  by  the  same  forces.  Working  our  way  backwards  from  the  j1*1  plane, 
we  locate  the  first  plane  whose  motion  differs  appreciably  from  that 
of  the  corresponding  plane  in  the  second  segment.  Thereby,  we  determine 
precisely  the  location  of  the  shock  front  at  any  time.  Furthermore, 

until  the  shock  reaches  the  2L./*1  plane,  it  is  necessary  to  solve  the 

equations  of  motion  for  only  particles  in  the  first  two  segments  since 
the  particles  in  all  the  remaining  segments  will  have  trajectories 
identical  to  the  corresponding  particles  in  the  second  segment.  When 

the  shock  front  nears  the  2L„  plane,  particles  in  the  third  segment 

are  included  in  the  computation.  The  shock  front  is  located  in  the 

same  manner  as  before  and  the  process  repeated  as  often  as  necessary  to 

complete  the  calculation.  It  is  therefore  necessary  to  monitor  at  most 

4l.„  planes  of  particles  ahead  of  the  shock  front  at  any  time,  the  last 

“ 2 
2L_  representing  equilibrium  conditions  ahead  of  the  shock  . 

3 . 2 Equations  of  Motion. 

The  preceding  subsection  contains  a qualitative  description  of  how 
the  lattice  is  subjected  to  shock  compression.  In  this  subsection  the 
equations  of  motion  to  be  solved  numerically  for  each  atom  in  the 
lattice  will  be  determined.  As  has  been  observed  before,  it  will  be 
most  convenient  to  use  nondimensional  variables  in  all  calculations. 

For  convenience,  Table  III  lists  the  normalizing  factor  for  various 
quantities  of  interest.  Hereafter,  reference  to  these  quantities  will 
imply  their  nondimensional  values. 


The  force  exerted  on  particle  (i,j,kj  by  the  remaining  particles 
in  the  lattice  can  be  determined  from  the  gradient  of  Eq.  (2.14)  with 
respect  to  If  we  normalize  the  force  by  2RD/ro  (see  Table  III) 

and  employ  nondimensional  units  for  the  time  and  for  R.jj.»  the  equa- 
tion of  motion  for  the  particle  becomes 


(3.5) 


2.  We  are  indebted  to  D.H.  Tsai  for  suggesting  this  method  of  calcu- 
lation to  us. 
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TABLE  III.  Normalizing  Factors. 


Quantity 

Symbol 

Normalizing  Factor 

Velocity  of  particle  (i,j,k) 

Vijk 

\A)/m 

Nondimensional  lattice  constant 

A 

o 

ro 

Position  of  particle  (i,j,k) 

Kijk 

a 

o 

Potential  energy  of  particle  (i,j,k) 

$ 

ijk 

D 

Kinetic  energy  of  particle  (i,j,k) 

T 

ijk 

D 

Density 

P 

po 

Force  exerted  on  particle  (i,j,k) 

?ijk 

2RD/r 

Time 

T 

\rtn/D  aQ 

Element  of  stress  tensor 

°ae 

AB/al 

Compression  velocity 

UP 

VD/m 

Temperature 

0 

D/k„ 

Dissociation  energy 

D 

Mass  of  particle 

m 

Equilibrium  spacing  for  isolated 
pair  of  atoms 

ro 

Lattice  constant 

a 

o 

Initial  density  ahead  of  shock 

po 

Here,  ?. ^ is  the  nondimensional  force  exerted  on  the  particle  given 
explicitly  by 

U(A0,V*tJ  •1)  -R<\Ajk4«.J 

lJk  ' J.n  [' 

1^...-^.  I<  R 

1 ijk  tmn'  c /•, 


X 


^ijk’^fcmn 

^ijk'^imJ 
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. 


- -«#*  . . 


and  each  dot  represents  differentiation  with  respect  to  the  dimension- 
less time,  t,  related  to  the  real  time,  t,  by 


(3.7) 


The  sum  in  Eq.  (3.6)  runs  only  over  values  of  t,m,  and  n such  that 

|$.  .. -ft.  | < R . 

1 ijk  tmn1  c 

For  purposes  of  computation  it  is  most  convenient  to  convert 
Eqs.  (3.5)  into  a set  of  first-order  equations.  This  transformation 
can  be  accomplished  by  noting  that  the  velocity  of  particle  (i,j,k)  is 
just  the  time  derivative  of  the  displacement.  Equation  (3.5)  becomes, 
therefore, 


Vijk  “2AoR  Fijk 


(3.8) 


ft. 

ijk 


Vijk 


We  have,  of  course,  twice  the  original  number  of  equations. 


For  the  equilibrium  lattice,  Eqs.  (3.8)  apply  to  every  particle 
in  the  lattice.  For  the  lattice  undergoing  shock  compression,  however, 
particles  in  the  first  plane  are  not  acted  on  by  any  net  force  in  the 
2 direction.  Therefore,  we  have 


( u 


u 


k = 1,2 


(3.9) 


where  is  the  constant  compression  velocity  and  the  subscript  2 


denotes  the  component  in  the  2 direction, 
in  the  lattice,  Eqs.  (3.8)  apply. 


For  the  remaining  particles 


In  order  to  solve  Eqs.  (3.8)  and  (3.9),  we  employed  a fourth-order 

Runge-Kutta  technique^.  To  check  the  accuracy  of  the  code  and  the  appro- 
priateness of  the  step  si2e  used  we  calculated  the  work  done  on  the 


B.  Carnahan,  H.A.  Luther,  and  J.O.  Wilkes,  Applied  Numerical 
Methods  (Wiley,  NY,  1969),  Chap.  6. 
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lattice  by  the  compression  force  and  compared  it  at  various  times  with 
the  increase  in  the  total  energy  of  the  lattice.  For  most  calculations, 
a step  size  of  the  order  ,05/R  was  found  to  be  acceptable. 


4.  THERMODYNAMIC  VARIABLES 


Once  Eqs.  (3.8)  have  been  solved  numerically  to  determine 
and  the  results  can  be  used  to  calculate  macroscopic  quantities 


of  interest  in  the  lattice.  These  thermodynamic  quantities  vary  with 
position  and  therefore  the  averages  from  which  they  are  calculated  are 
taken  only  over  finite  regions  of  the  crystal.  Hereafter,  we  will  re- 
fer to  such  regions  as  "thermodynamic  regions".  We  will  assume  that 
the  regions  are  sufficiently  small  that  the  variable  of  interest  varies 
negligibly  throughout  the  region,  and  yet  sufficiently  large  that  the 
averages  are  meaningful  and  boundary  effects  negligible.  For  all 
thermodynamic  variables,  standard  statistical-mechanical  definitions 
will  be  employed  and  ensemble  averages  will  be  replaced  by  spatial 
averages.  In  the  following  discussion,  we  will  assume  that  the  ther- 
modynamic region  over  which  the  average  is  taken  begins  just  before 
some  odd-numbered  plane  (normal  to  the  z axis)  and  ends  just  before  a 

later  odd-numbered  plane.  Respective  values  of  k are  k . and  k 

r r min  max 

In  the  y and  z directions  the  region  is  defined  by  the  same  bounds  as 
the  primary  lattice.  Therefore,  in  calculating  averages  in  a thermo- 
dynamic region  we  include  those  particles  whose  indices  lie  within 


1 < i < L 


1 < j < L 


k . < k < k 

min  max 


(4.1) 


The  total  number  of  atoms  which  satisfy  conditions  (4.1)  is  given  by 


nD  = (k, 


max 


- k 


nun 


1) 


L L . 
x y 


(4.2) 


4.1  Average  Velocity. 


The  shock  will  induce  a flow  in  the  direction  of  its  propagation 
and,  thus,  it  is  necessary  to  calculate  the  average  velocity  as  a 
function  of  position  in  the  lattice.  For  the  thermodynamic  region 
specified  by  Eq.  (4.1)  the  appropriate  expression  is 


< $ > 


1 


^max 

I 

k»k_ 


(4.3) 
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where  £ is  a short-hand  representation  indicating  that  the  sun  extends 
R 

only  over  the  atoms  contained  in  the  region.  One  would  expect  that  only 
<VT>  would  be  nonzero. 

4. 2 Density. 

Prior  to  the  time  the  lattice  is  compressed  by  the  shock  wave,  its 

3 

density  is  a constant  and  equal  to  4/a\  This  value  can  be  obtained 

from  Figure  1 if  we  note  that  one-eighth  of  each  particle  at  the  corners 
of  the  cube  and  one-half  of  each  particle  at  the  faces  of  the  cube  are 
contained  within  the  volume  of  the  cube.  Compression  by  the  shock  does 
not  affect  the  average  separation  between  atoms  in  the  x and  y direc- 
tions, but  obviously  does  in  the  z direction.  Thus,  the  final  density, 
normalized  by  the  original  density,  can  be  obtained  by  calculating  the 
final  separation  along  the  z axis  of  planes  associated  with  k^^  and 

K and  dividing  by  their  equilibrium  separation.  In  the  dynamic  case, 

the  z coordinate  of  a plane  will  be  defined  as  the  mean  z component  of 
the  atoms  within  it.  Consequently,  from  F.q.  (2.3)  the  initial  separa- 
tion of  the  planes  defining  the  region  in  question  is 


1’k  -ll  fk  . -1  1 ) 
-I--  - m-y-  | , 


(4.4) 


and  the  final  separation  is 


*Vrr  \ X’k.j., 

x y i = l j =1  [ 


'i.j.k  -1 
J max 


(4.5) 


‘i.j.k  . “i.j.k  . +1 
nun  min 


The  final  density,  therefore,  is 


(4.6) 


4.3  Temperature. 


The  usual  definition  of  temperature,  normalized  by  the  factor 
D/kg , is  given  by  the  expression 

'■A  l 


(4.7) 


It  is  also  convenient  to  define  a "component"  of  the  temperature  in 
each  Cartesian  direction.  A reasonable  definition  is 


9 - -r—  1 |(V.  ..)  - < V >T 

“ 3nR  R \ ljk  a “ J 


(4.8) 


4.4  Stress. 


The  ati  component  of  the  stress  tensor  in  a solid  is  defined  as  the 

ath  component  of  the  force  acting  per  unit  area  on  a plane  normal  to 
the  6 direction.  Such  forces  arise  both  from  the  interactions  acting 
across  the  plane  as  well  as  from  transfer  of  momentum  due  to  fluctua- 
tions in  the  velocity  of  the  atoms  about  the  mean  velocity.  For  dilute 
gases,  the  latter  term  is  dominant  and  represents  the  pressure  tensor. 

For  solids  or  dense  gases,  however,  the  potential  contribution  must  be 
accounted  for. 

4 

Irving  and  Kirkwood  have  derived  an  expression  for  the  ensemble- 
averaged  pressure  tensor  for  a system  in  which  the  potential  energy  is 
not  negligible.  For  purposes  of  computer-molecular-dynamic  calculations, 
we  assume  that  the  ensemble  average  can  be  approximated  by  a spatial 
average  over  the  appropriate  thermodynamic  region.  Provided,  then,  that 
the  range  of  the  interatomic  forces  is  much  less  than  the  dimensions  of 
the  region  in  question,  the  kinetic  and  potential  contributions  to  the 
pressure  tensor  become  [see  Reference  4,  Eqs.  (5.13)  and  (5.14)] 


— y tv...  -<?>)tf..„  -<v>) 

n„  U 1 ilk  ilk  ’ 


(4.9) 


’V  * "R  R ti.n  «■" 


(4.10) 


where  ^jjj'tmn  *s  torce  exerte<i  on  particle  (i,j,k)  by  particle 

(t,m,n).  In  transcribing  the  equations  from  Reference  4,  we  have  em- 
ployed our  nondimens ional  units  and  normalized  the  stress  by  the  fac- 
tor 4D/a^  (see  Table  III).  The  total  stress,  of  course,  is  given  by 


°k  * °V 


(4.11) 


71  J.H.  Irving  and  J.G.  Kirkwood,  "The  Statistical  Mechanical  Theory 
of  Transport  Processes.  IV.  The  Equations  of  Hydrodynamics",  J. 
Chem.  Phys.  18,  817  (1950). 


24 


ArPKNDIX 


In  this  appendix  we  briefly  describe  the  computer  program  which 
performs  the  calculations  described  previously.  The  description  is 
followed  by  a list  of  the  major  symbols  used  in  the  code,  their  defini- 
tions, and  the  corresponding  symbol  used  in  the  text.  In  Section  A. 2 
we  describe  the  function  of  each  of  the  sixteen  subprograms  in  the 
program.  Finally,  a complete  listing  of  the  code  is  presented  in 
Sect  ion  A. 3. 

The  program  consists  of  a main  routine  and  sixteen  subprograms  which 
will  be  discussed  in  the  following  section.  It  is  capable  of  calculating 
a set  of  initial  conditions  for  the  primary  lattice  prior  to  shock  com- 
pression and  of  performing  the  shock-wave  calculation.  At  the  end  of 
each  run  the  velocities  and  positions  of  each  atom  as  well  as  the  time 
are  written  on  a tape.  Thus,  by  using  these  values  as  input  data,  a 
succeeding  calculation  can  begin  where  a previous  one  was  terminated. 

The  execution  time  for  the  program  is  approximately  7 X 10  ^ seconds 
per  particle  per  time  step  on  a CDC  7b00. 

The  basic  function  of  the  program  is  to  solve  numerically  the 
classical  equations  of  motion  for  every  atom  in  the  lattice  and,  from 
their  solution,  to  calculate  a number  of  quantities  of  interest.  Includ- 
ed in  the  possible  output  are  the  mean  displacements  and  velocities  of 
planes  of  atoms  normal  to  the  z axis  at  a particular  time;  the  velocity- 
time trajectories  of  individual  planes;  the  velocity  distribution 
function  and  various  thermodynamic  quantities  calculated  in  macroscopi- 
cally  small,  microscopically  large  regions  of  the  lattice  (referred  to 
as  thermodynamic  regions);  and  the  velocity  and  displacement  of  each 
atom  in  the  lattice  at  a given  time. 

A number  of  input  data  must  be  supplied  to  the  main  routine.  These 
data  are  listed  and  described  as  follows: 

LX  (statement  10)  - one-half  the  total  number  of  planes,  normal  to  the 
x axis,  of  the  primary  lattice 

LY  (statement  10)  - one-half  the  total  number  of  planes,  noi'mal  to  the 
y axis,  of  the  primary  lattice 

L2  (statement  10)  - one-half  the  total  number  of  planes,  normal  to  the 
z axis,  whose  atomic  equations  of  motion  must  be  solved  in  the 
current  calculation 

NN1AX  (.statement  14)  - the  final  set  of  equilibrium  nearest  neighbors  to 
a particle  which  are  scanned  to  determine  whether  they  lie  within 
the  critical  radius  of  potential  interaction. 

NMIN  (statement  14)  - the  last  set  of  equilibrium  nearest  neighbors  to 
a particle  which  are  assumed  to  lie  within  the  critical  radius  of 
potential  interaction.  Therefore,  if  NMIN=2  and  NMAX=3,  for  example, 
all  particles  which  are  first-  or  second-nearest  neighbors  to  a 
particular  particle  when  the  lattice  i->  in  equilibrium  are  assumed 
to  lie  within  the  critical  radius  in  the  dynamic  lattice;  those 


\ 
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which  are  fourth-  or  higher-nearest  neighbors  in  equilibrium  are 
assumed  to  be  outside  the  critical  radius;  those  which  are  third- 
nearest  neighbors  are  tested  to  determine  whether  they  are  within 
the  critical  radius. 

RCR1T  (statement  14)  - the  radius  of  the  imaginary  sphere  drawn  about 
each  particle  whose  equations  of  motion  are  solved.  Atoms  lying 
within  the  sphere  are  assumed  to  exert  a force  on  the  given  atom; 
those  outside  do  not. 

H(statement  15)  - the  time  step  employed  in  the  Runge-Kutta  scheme  for 
numerically  solving  the  equations. 

TAUMAX  (statement  15)  - the  time  at  which  present  calculation  is  to  stop. 

IPRINT  (statement  16)  - the  frequency  of  printout  of  all  output  data 
except  velocity-time  trajectories  of  specific  planes.  That  is, 

IPRINT  time  steps  are  executed  between  printouts. 

R (statement  17)  - anharmonicity  factor  in  the  Morse  potential 

RE  (statement  17)  - lattice  constant,  o.Q,  normalized  by  single-pair, 
equilibrium  separation  constant,  r . 

GAMMA  (statement  18)  - ratio  of  dissociation  energy,  D,  in  Morse  potential 
and  thermal  energy,  kfiT,  of  lattice. 

ICALC  (statement  21)  - inc-x  which  determines  which  of  three  calcula- 
tions is  to  be  performed  by  the  code.  For  ICALOl,  the  program 
will  calculate  random  initial  conditions  for  the  atoms  in  the  pri- 
mary lattice  prior  to  shock  compression;  for  ICAL02,  the  program 
reads  this  initial  data  from  TAPE  7 and  performs  the  first  shock- 
wave  calculation;  for  ICAL03,  the  program  reads  from  TAPE  7 the 
positions  and  velocities  of  all  atoms  in  the  lattice  as  well  as 
the  time  at  the  end  of  the  last  shock-wave  calculation,  and  extends 
the  shock-wave  calculation  to  TAUMAX. 

K21  (statement  22)  - the  number  of  planes,  normal  to  the  z axis,  which 

make  up  a thermodynamic  region.  Thus  if  K21=8,  for  instance,  thermo- 
dynamic properties  of  the  lattice  will  be  calculated  for  planes 
1-8,  9-16,  etc. 

The  remaining  parameters  are  input  data  only  for  ICALC* 1: 

LZLAST  (statement  56)  - one-half  the  number  of  planes,  normal  to  the 
z axis,  contained  in  the  previous  calculation.  This  will  differ 
from  LZ  if  one  must  increase  the  size  of  the  lattice  to  perform 
the  present  shock-wave  calculation  (see  discussion  Section  3.1). 

LZSEG  (statement  56)  - one-half  the  number  of  planes,  normal  to  the 
z axis,  contained  in  the  primary  lattice. 

NSEG  (statement  56)  - the  number  of  segments  equal  in  size  to  the  pri- 
mary lattice  which  must  be  added  to  the  lattice  of  the  previous 
calculation  to  carry  out  the  current  calculation,  i.e.,  LZ«LZLAST 
+ NSEG*LZSEG. 
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UP  (statement  57)  - compression  velocity  (normalized  by  vb/"m  ). 

NPMAX  (statement  58)  - the  total  number  of  planes  whose  velocity-time 
trajectories  are  to  be  printed. 

KPRINT  (statement  58)  - frequency  of  printout  for  velocity-time  tra- 
jectories of  particular  planes  . 

NUMPL  (l)  (statement  61)  - number  of  particular  plane  whose  velocity-time 
trajectory  is  to  be  printed  out  (1  < I < NPMAX). 

DELT(I)  (statement  61)  - the  time  interval,  beginning  at  time  TAUI(I). 
for  which  the  velocity-time  trajectory  of  plane  number  NUMPL(I)  is 
printed. 

TAUI(I)  (statement  61)  - time  at  which  printout  of  trajectory  of  plane 
NUMPL(I)  is  begun. 


A.l.  List  of  Variables. 

This  section  contains  in  tabular  form,  a list  of  important  symbols 
employed  in  the  code,  their  definitions,  and  the  corresponding  symbol 
used  in  the  text. 


TABLE  A.l. 

Variables  in  Computer  Program. 

Symbol 

Symbol 

in  Code 

Definition 

in  Text 

DELT(I) 

Time  interval  for  trajectory  printout  of  plane 
NUMPL (I) 

DENS 

Particle  density  in  thermodynamic  region  of 
lattice 

P 

DV 

Velocity  interval  in  calculation  of  distribution 
function 

DXAV 

Average  x component  of  displacement  of  particles 
in  a plane 

DYAV 

Average  y component  of  displacement  of  particles 
in  a plane 

DZAV 

Average  z component  of  displacement  of  particles 
in  a plane 

ECHECK 

Sum  of  initial  energy  in  lattice  and  total  work  i 

done 

E1NIT 

Initial  energy  in  primary  lattice 

EP 

Potential  energy  of  a single  particle 

♦ . .. 
ijk 

TABLE  A. 1 . (Continued) 

Symbol 
in  Code 

Definition 

Symbol 
in  Text 

EPAV 

Average  potential  energy  of  single  particle 
in  thermodynamic  region  of  lattice 

EPOT 

Total  potential  energy  in  thermodynamic  region 
of  lattice 

ETHERM 

Thermal  (kinotie)  energy  in  thermodynamic  region 
of  lattice 

ETOT 

Total  energy  in  lattice 

F.TRANS 

Translational  (kinetic)  energy  in  thermodynamic 
region  of  lattice 

FX 

x component  of  force  exerted  on  a particle 

<PiJk>x 

FY 

y component  of  force  exerted  on  a particle 

(Fi.tkV 

F2 

z component  of  force  exerted  on  a particle 

<Fijk>* 

liAMMA 

Ratio  of  dissociation  to  thermal  energy  in 
lattice 

Y 

H 

Time  step 

I CALC 

Index  to  determine  which  calculation  code  is 
to  perform 

I PR I NT 

Index  to  determine  printout  frequency  for  output 

K21 

Number  of  planes  (normal  to  z axis)  contained 
within  a thermodynamic  region 

LX 

One-half  the  number  of  planes  of  atoms  in  the 
lattice  normal  to  the  x axis 

L 

X 

LY 

One-half  the  number  of  planes  of  atoms  normal 
to  the  y axis 

L 

y 

LZ 

One-half  the  number  of  planes  of  atoms  normal 
to  the  z axis 

LZLAST 

One-half  the  number  of  planes  of  atoms,  normal 
to  the  z axis,  in  the  previous  calculation 

LZSEG 

One-half  the  number  of  planes  of  atoms,  normal 
to  the  z axis,  in  the  primary  lattice 

L 

Z 

NK1K2 

Total  number  of  atoms  in  thermodynamic  region 

nR 

NMAX 

Largest  set  of  equilibrium  nearest  neighbors 
scanned 
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TABLE  A. 1.  (Continued) 


Syrabo l 
in  Code 

NM1N 

NP 

N PLANT. 
NPMAX 

NSEC 

NTOT 

NIIMPL(I) 

PKXX 

PKYY 

pkz: 

PVXX 

PVYY 

pvzz 

PXX 

PYY 

PZZ 

R 

RE 

T 

TAll 


Definition 


Symbo  1 
in  Text 


Largest  sot  of  equilibrium  nearest  neighbors 
assumed  to  bo  within  RCRIT 

Number  of  particles  per  plane  normal  to  the 
z axis  in  primary  lattice 

Number  of  plane  normal  to  the  z axis 

Total  number  of  planes  whose  trajectories 
are  printed 

Number  of  segments,  equal  in  size  to  the 
primary  lattice,  added  to  lattice  of  previous 
calculation,  i.e.  LZ»LZLAST*NSEC*LZSEG 

Total  number  of  particles  in  lattice 

One  of  NPMAX  planes  whose  trajectories  are 
printed 

xx  component  of  kinetic  contribution  to 
pressure  tensor 

yy  component  of  kinetic,  contribution  to 
pressure  tensor 

component  of  kinetic  contribution  to 
pressure  tensor 

xx  component  of  potential  contribution 
to  pressure  tensor 

yy  component  of  potential  contribution 
to  pressure  tensor 

zz  component  of  potential  contribution 
to  pressure  tensor 

xx  component  of  pressure  tensor 
yy  component  of  pressure  tensor 
zz  component  of  pressure  tensor 
Anharmonicity  factor  in  Morse  potential 
Lattice  constant,  normalized  by  single- 
pair equilibrium  separation,  r^,  in  Morse 
potential 

Temperature  in  thermodynamic  region  of  lattice 
Time 


n 


P« 


(a.  ) 

k xx 


(Vyy 

<v« 


(o  ) 

V XX 


(o  ) 
v yv 


(Vzz 


0 

XX 


0 

T 


TABLE  A.l, 

(Continued) 

Symbo 1 

Symbol 

in  Code 

Definition 

in  Text 

TAUI(I) 

Time  at  which  printout  of  trajectory  of  plane 
NUMPL(I)  is  begun 

TAUMAX 

Maximum  time  for  which  calculation  is  carried  out 

TX 

x component  of  temperature  in  thermodynamic  region 

6 

X 

TY 

y component  of  temperature  in  thermodynamic  region 

0 

y 

TZ 

z component  of  temperature  in  thermodynamic  region 

0 

z 

UP 

Compression  velocity 

U 

p 

VXd.J.K) 

x component  of  velocity  of  particle  (i,j,k) 

fVijk>x 

VY(I  ,J,K) 

y component  of  velocity  of  particle  (i,j,k) 

(vi)kV 

VZ(I.J.K) 

z component  of  velocity  of  particle  (i,j,k) 

(Vijk>z 

VXAV 

Average  x component  of  velocity  of  particles 
in  a plane  normal  to  the  z axis 

VYAV 

Average  y component  of  velocity  of  particles 
in  a plane  normal  to  the  z axis 

VZAV 

Average  z component  of  velocity  of  particles 
in  a plane  normal  to  the  z axis 

X(I.J.K) 

x component  of  position  of  particle  (i,j,k) 

xijk 

YO.J.K) 

y component,  of  position  of  particle  (i,j,k) 

Yijk 

Z(I.J.K) 

z component  of  position  of  particle  (i,j,k) 

Zijk 

A . 2 . Description  of  Subprograms. 

This  section  contains  a list  of  the  sixteen  subprograms  with  a 
basic  description  of  the  function  of  each  of  the  routines.  The  list 
(alphabetical)  is  as  follows: 

CALDIST  - The  subroutine  calculates  the  components  of  the  distance 
between  two  particles  whose  potential  energy  or  force  of  inter- 
action is  being  computed.  The  purpose  is  to  determine  whether 
one  particle  lies  within  the  critical  radius  of  the  other. 

CONVERT  - The  subroutine  converts  integers  (i,m,n),  identifying  a 

particle  outside  the  primary  lattice,  to  the  corresponding  atom 
inside  the  primary  lattice  identified  by  indices  (LEQ,MEQ,NEQ) . 

DISTFN  - The  subroutine  calculates  the  velocity  distribution  function 
in  thermodynamic  regions  of  the  lattice.  More  precisely,  it 
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determines  the  fractional  number  of  atoms  whose  velocity  at  any 
time  has  components  lying  between  VX  and  VX+DV,  VY  and  VY+DV, 
and  V Z and  VZ+DV. 

FORCE  - The  subroutine  calculates  the  components  of  the  force  exerted 
on  one  par*ide  by  another,  the  distance  between  the  atom  being 
provided  by  CALDIST. 

NABORS  - The  subroutine  calculates  the  differences  I-L,  .J-M,  and  K-N 
for  any  particle  (I,J,K)  and  its  neighbor  (L,M,N).  The  values 
are  stored  in  an  array  for  future  use  in  the  program. 

NRAN31  - The  subroutine  generates  random  numbers,  consistent  with  a 
Gaussian  distribution,  for  assigning  velocities  to  atoms  of  the 
primary  lattice. 

OUTPUT  - The  subroutine  prints  the  single  - particle  displacements  and 
velocities  for  every  atom  whose  equations  of  motion  are  solved. 
Also  printed  are  the  total  sums  of  x,y,  and  z components  of 
velocities . 

PLANEAV  - The  subroutine  averages  the  x,y,  and  z components  of  dis- 
placements and  velocity  for  every  atom  in  a single  plane  and 
prints  the  results.  Also  printed  are  the  total  sums  of  the 
components  for  all  atoms  in  the  lattice. 

POTFOR  - The  subroutine  calculates  the  potential  energy  of  a single 
particle  and  the  contribution  to  the  potential  part  of  the 
stress  tensor  from  that  particle. 

SAVE  - The  subroutine  saves  current  values  of  position  for  each  atom 
in  the  vector  (XS(I,J,K),  YS(I,J,K),  ZS(I,J,K)). 

SCAN  - The  subroutine  determines  which  neighbors  between  NMIN  and 
NMAX  lie  within  the  critical  radius  of  a given  particle  and 
calls  the  subroutines  which  calculate  the  potential  of  the 
particle  as  well  as  the  force  exerted  on  it.  For  the  first 
iteration  of  the  Runge-Kutta  scheme,  the  nearest  neighbors 
are  calculated  and  are  written  on  TAPE  8;  for  subsequent 
iterations  atoms  within  the  critical  radius  are  assumed  not 
to  change. 

SOLVE  - The  subroutine  uses  a fourth-order  Runge-Kutta  scheme  to 
numerically  solve  the  equations  of  motion  for  each  atom  in 
the  lattice.  It  also  calculates  the  work  done  on  the  lattice 
by  the  external  driving  force  in  a particular  time  interval. 

SORT  - The  subroutine  (used  in  conjunction  with  DISTFNj  counts  the 
number  of  atoms  lying  in  each  velocity  interval. 

START  - The  subroutine  assigns  atoms  to  their  initial  equilibrium 

positions  and  generates  their  initial  velocities  consistent  with 
a Maxwellian  distribution. 


THERMO  - The  subroutine  calculates  the  temperature,  density,  pressure 
tensor  and  mean  velocity,  averaged  over  thermodynamic  regions  of 
the  lattice. 

URAN31  - This  function  subprogram  is  called  from  NRAN31  and  used  to 
generate  random  numbers  for  initial  conditions  in  the  primary 
lattice.  r ’ 


A.  3.  Listing  of  Computer  Program. 

The  final  section  of  this  appendix  contains  a complete  listing  of 
the  main  program  and  each  of  the  sixteen  subprograms. 


00011 

0002: 

0003! 

0004  i 

0003! 

0004! 

0007! 

0000! 

0003! 

0010! 

0011! 

0012! 

0013! 

0014  ! 

0013! 

0014! 

0017! 

0010! 

0019! 

0020! 

0021! 

0022! 

0023! 

0024! 

0023! 

0024! 

0027! 

0020! 

0029! 

0030! 

00311 

00321 

0033!' 

00341 

0033! 

0034! 

0037! 

00301 

00391 

0040! 

0041! 

0042! 

0043! 

0044! 

0043! 

00441 

0047! 

0040! 

0049! 

0030! 

0031! 

0032! 

0033! 

0034! 

0033! 


PROGRAN  BATTCHI INPUT. OUTPUT. T APES* INPUT. TApC4»0UTPUT«T APCT. TAPES  I 
CONNON/tfCL/VXI04,04«300).VT(04.04,300) ,V2|04,04,300) 

C OMMON/POS I TS/XS I 04, 04. 300 > .TS( 04, 04. 300) .23104.04.3001 

CONNON/PUSI T/X I 04, 04. 300). 7(04, 04. 3«*0 1 ,2(04.04,3001 

CONNON/NNBORS/IX (4.24,2 1 , It (4.24.2) ,12(4,24.2 >,NUNt4> 

CONNON/SEARCH/NHAX.NHIN,NN,RCAlT.ACAIT2 

COMNON/S 1 2E/LX , LT ,Li*  NP .NT0T.U24.L2 SC 64 

COMMON/M t SC/GAMMA, TAU.H.FX, FT .F2.R.RC. CP. 1CACC.CCHCCR 

01NCNSI0N  NUMPL (201,11(20), TAUK20I , OCLT (201.TAUI120) 

RCADIS,  1001  LX, (.7 ,1.2 

C LX.LT.ANO  L2  ARC  OINCNSIONS  THIS  CALCULATION.  LX.L7.  AND  L2LAST  ARC 
C OINCNSIONS  OF  LAST  CALCULATION.  LX.LT , AND  L2SC6  ARC  OINCNSIONS  FOR  A 
C UNIT  SC6NLNT. 

RCADIS, 101)  NNAX.NNIN.RCRIT 
RCADIS. 102)  H.TAUNAX 
RCAO(3« 103 1 I PRINT 
RCAOIS. 104)  R.RC 
RCAOIS, 109)  GANNA 

C ICALC*1«2.0R3  FOR  INlTIAL>CONOlTION  CALCULATION. INITIAL  SMOCK-WAVC 
C CALCULATION. OR  SUBSCOUCNT  SHOCK-WAVE  CALCULATIONS, RCSPCCTIVCLt. 
RCAOIS, 104)  ICALC 
RCAOIS. 1071  K21 

C K21  IS  THC  NUNBER  OF  PLANCS  OVCR  WHICH  THC  DISTRIBUTION  FUNCTION  IS 
C CALCULATED. 

NINT*2*LZ/K21 

NP«2*LX*LT 

L24*4*L2 

NT0T*NP*L24/2 

RCR1T2*RCRIT**2 

IF (ICALC  .CO.  II  L2SCG*L2 

IF (ICALC  .CO.  1)  L2SE64*4«c2SCG 

WRITE  1* , 10#  ) LX.LT.L2 

WRITCI4.109I  NNAX.NNIN.RCRIT 

WR1TCI4.I10)  H.TAUNAX 

WRITE (4.111)  R.RC 

WRITCI4.112I  GANNA 

CALL  NABORS 

IF (ICALC  «NC , II  GO  TO  10 
CALL  START 

C CALCULATE  2CRO- POINT  POTENTIAL  CNCR67 . 

CALL  SAVE 

CALL  SC AN (LX.LT, L24> 2,1) 

REWIND  • 

CPT0TwEP4FL0ATINT0T! 

C CALCULATE  INITIAL  KINETIC  CNCR6T  IN  LATTICC. 

CKTOTwO.O 
00  4 1*1, LX 
00  4 J*1 «LT 
DO  4 K*1 «L24 

4 CKTOT»CKTOT ♦0,3*(VXII.J.K)**2*VT(I,J,K)**2»V2(I,J.KI**2I 
CCHCCR*CKTOT*CPTOT 
WRITCI4.113)  CPTOT.CKTOT.CCHCCK 
TAU*0.0 
GO  TO  43 

C BCG  IN  SCOUCNCC  FOR  PCRFORNlNG  INITIAL  SHOCK-WAyC  CALCULATION. 


oo3*: 
0097 : 

oooo : 

0099 1 
00401 
00411 
00421 
0049! 
00401 
0049: 
0044: 
0047: 

0040 : 

0049: 

0070: 
0071 : 

0072: 

0079: 

007*: 

0079: 

0074: 

0077: 

oo;o: 

0079: 

oooo: 
oooi: 
0002: 
ooo9: 
ooo*: 
0009: 
oooo: 
0007: 
oooo: 
0009: 
0090 : 

00911 

0092: 

00991 

009*: 

00991 

00941 

00971 

00901 

0099.* 

01001 

0101*. 

01021 

0109: 

010*1 

0109) 

0104: 

01071 

01001 

01091 

01101 


10  RCAU19.il* I L2LA9T.L2SCG.N3C0 
ACAD(9.119)  UP 
RCA0I9.120I  NPNAX.KPAINT 
IF1NPNAX  .CO.  0>  00  70  19 
00  12  LLal.NPMAX 

ACA0<9.121INUNPL(LD.0CLT1LL>.TAUI1LD 
IIILLtaO 
12  TAUliLLIaO.O 
19  WRITCtO.llOl  UP 

WRITC14.119)  L2SEG.L2LAST 
L2L AS*«LZL AST •* 

L2SC6*aLZSCG«* 

00  19  I«1.LX 
DO  19  J«1*LY 
00  19  K*1.L2LAS« 

19  ACA0I7)  XII. J.K)  .Y(I.J.K),ZlI.J,K>,VX<l,J,KI.VYtI.J,KI,VZ(I,J,K) 
IFIICALC  .CO.  91  GO  TO  *0 

REWIND  7 
TAUaO.O 

00  29  llal iNSCG 
00  20  I»1.LX 
DO  20  J*1 «LY 
00  20  M1.L2SCG* 

KPaIULZSCG*»K 

XlI.J.KPlaXlI.J.K) 

VII.J.KPIaril.J.K) 

Z 1 1 » J.KP)aZ< I * J »K ) +0 ,9»FL0aT ( lKP*l ) /2I *«9*FL0AT1 (K*l 1/21 

VXU.J.KPtaVXtl.J.KI 

VTtl.J.KPIaVYlI.J.KI 

20  VZlI.J.KPIaVZlI.J.K) 

29  CONTINUE 

00  30  lal.LX 
DO  90  jal.LY 
VZU'Jdl«UP 
90  V2 ( I ♦ J.2 )*UP 

C CALCULATE  INITIAL  CNEROY  IN  LATTICE. 

CKTOTaO.O 
CPTOTsO.O 
00  99  lal.LX 
DO  99  jaltLY 
00  99  Kal.LZ* 

CXTOT»CKTOT*0 . 3* I VX II«J.K)«*2*VYt I , J.K I **2«V2( I , J.K) **2> 

CALL  SAVE 

CALL  SCANI I .J.K .2  < 1 > 

99  CPTOT»CPTOT*EP 
REWIND  0 

ECMCCKaCPTOUCKTOT 
WRITC14.U7)  TAU.CCHECK 
CALL  PLANCAV 
00  TO  49 

*0  RCA0I7I  TAU.CCHECK, CT0T1 
RCN1N0  T 

IFlNSCO  .CO.  0>  00  TO  40 
00  90  Ilal.NSCO 
00  *9  lal.LX 
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I 


4 


i 


11 

12 

13 

14 

15 
IS 
17 
IS 

19 

20 
21 
22 
23 
29 
25 
20 
27 
20 

29 

30 

31 

32 

33 

34 
33 
30 
37 
30 

39 

40 

41 

42 

43 

44 
43 
40 
47 
,40 

49 

30 

31 

32 

33 

34 
33 

50 
37 
30 
59 
00 
01 
02 
03 
04 
03 


00  43  Jal,  V 
00  43  K*1 . LZSE64 
KPbL2LAS4» 111*11 *L2SC04*K 
KPP»LZLAS4-L2SCG4*K 
Kl 1 » J.KP ) aXI I • J.KPP) 

T 1 I .J.KP I at  1 1 .J.KPP I 

2 1 1 » J • KP I »Z 1 1 • J • KPP | ♦ 0.S*FL0AY< IKP-1I/2I-0.94FLOOTHKPP-1I/2* 
VXII.J.KPIaVXII .J.KPP I 
VYII.  J.KPIWYI  I.  J.KPP) 

43  V2II.J.I«*IaVZII.J.KPPI 
30  CONTINUE 

c calculate  energy  in  lattice. 

EKTOTaC.O 
CPTOTaO.O 
DO  33  lal.LX 
DO  S3  J«1.LY 
00  33  Kal,L24 

EKTOTaEKTOT ♦0.3*tVXII.J»K)«*2*VYl 1 «J.K) a*2*VZt 1 .J.K I *at> 

CALL  SAVE 

CALL  SCANII.J.K.2.11 
35  EPTOTaEPTOT »EP 
REWIND  0 

ECHECRat CHECK aEPTOTaEKTOT'E TOT 1 
00  CONTINUE 

WRXTEI0.117I  TAU.ECHECK 
CALL  PLANEAV 
03  M«0 
70  *■*♦! 

1FITAU  .GE.  TAUNAX.H/2.01  60  TO  90 
CALL  SOLVC 

IFINPNAX  .EG.  01  60  TO  79 

DO  74  LLal.NPNAX 

IF ( T AU.LT .T AUI ILL  > >60  TO  74 

1 1 (LLI a 1 1 (LL  I ♦ 1 

TAU1 ILL I*TAU1 ILL  I AH 

IFITAUKLLI.CT.OELTILLII  00  to  74 

IFIIIILLI.NE.KPRINTI  GO  TO  74 

KKaNUHPLILLI 

KK2*2*KK 

KK2Hla2»KK-l 

VZAVaQ.O 

VXAVao . 0 

VTAVao.O 

00  72  lal.LX 

00  72  ja  ,LT 

VXAV«VXAV«VXI I • J.KK2 l«VX< I , J.KK2N1 I 
VYAVavYAV*VY<I.J.KK2)»VY<l . J.KK2M1 I 
72  VZAVBV2AV»VZtI.J.KK2l*V2U. J.KK2N1) 

V2 AVa V2 AV/FLOAT INP I 
VXAVavXAV/FLOATINP) 

VYAV»VT AV/FLOAT INP I 
VAV*SuRTIVXAV»»2*VYAV»*2*VZAV««2l 
WRITE 1 0.122 1 KK.V2AV.TAU.VXAV.VYAV.VAV 
IIILL>aO 
74  CONTINUE 


I < 


t 

I 


75  IF IM  ,N£.  IPHINTI  60  TO  TO 
CTOTaO.O 
MO 

UNITCI6. 117 1 TAU.CCHCCK 
IFUCALC  .NC.  1>  60  TO  00 
CALL  OUTPUT 
CALL  OUTPUT 
00  CALL  PLOW  A V 
00  05  l»l «NINT 
Kla!I.l)*K21«l 
K3al*K21 

CALL  0ISTFNIK1.K2.0.1) 

05  CALL  |HCNN0!K1,K2.CT0TI 
UNITE !6. 110)  CTOT 
00  07  Jal.LX 
00  07  jal.LT 
DO  07  Ml. LI* 

07  UNITE! T)  X(IfJ.K).T!l.J.K),Z!l.J.K).VX!l.J,K).VT(I,J.K)iVl!I<J,K) 
UNITE! 7)  TAU.CCMECH.CTOT 
NEU1N0  7 
60  TO  70 
00  STOP 
00  F0NNATIS1S) 

01  FONNAT II3.I3.F7.3) 

02  FONNAT I2F0.5 I 
03  FONNAT! 1*1 
00  FONNAT UFO, 31 
05  FONNAT !F6,S) 

06  FONNAT! 13 1 
07  FONNAT! IS) 

00  FONNAT !2X.NLXaO, l3.2X.OLVaO. 13.2X.OLlaO.I3) 

00  FONNAT ! IX.ONNAXaO, II . 3X.0NhIN*0, 12.SX.0NCNlTaN.E10.6l 

10  FONNAT I lX.8Mai,tIo.6»3X.OTAUNAxaO»El0.6) 

11  FONNAT I lX.0Na0,E10.6.3X.0NEa0«Cl*.6) 

12  FONNAT ! lX.06ANNAao.C10.6l 

13  FONNATI/IX.OZCNO'POINT  P0TCNTlALai.Cl*.4.3x.0INlTIAL  KINETIC  CNCN6 
2 T*O.C10.6.3X.OCCNCCM#.£10.6) 

10  FONNAT ! 3 1 0 I 

15  FONNATIFO.S) 

16  FONNAT ! lX.OUPaO.ElO.6) 

17  F0NNAT!lH1.0TAUaa.Cl*.6.3X.0CCHCCKaa,C10.6) 

10  FONNAT !lX.0CT0Ta0.Cl*.6l 

10  FONNAT ! IX .6LI3C6*0. 1 3 . SX.OLZLASTaO. IS  I 
20  F0NNATT2I3) 

22  FQNNAT!1X.NPL  NIJMO. 13 .5X,oy£AV«N.  C10.4.3x.0TAUa0,Cl0.4.S&« 
20VXAVa0.Cl0.6i3l.0VTAva0.Cl0»k.2Xi0vAVa0.Cl*.6) 

121  FONNAT! IS. 2F0. 51 
CNO 


3b 


Of  O#  Of  «V 


SUBROUTINE  CALDISTU'N>N.LCS.NCS.NC0<I.J.K) 

C0NN0N/STRCSS/DISTX.0ISTY.0ISTZ.013T.0IST2,C1.PVXX,PVYY.PVZZ 
CONMON/POSITS/X(0N.0N,3OO>,V(0N.0N,s00>.Z<0H.0N.300» 
IFU.CO.LCOIGO  TO  209 
IPPMUOtN.Nl 

IPsMOO  t NEO « H I /2*H00 • H* IPP  *H 1/2 

0ISTX«X<I.J.K)-X(LC0*N£0,NE0»*.9*FL0ATI2»LC8-2*LMPI 
60  TO  206 

209  OISTX»XU.  J.K»-XU.E8.NE0,NC0» 

206  IF (N.CB.MCOIGO  TO  219 

IP>MOU<  NEU-1 «2 I * 1 A0S  t N00(N-1 >2 II 

OISTYcVI I. J.KI-Y ILCe.NC0,NC0l».9*FL0AT(2*MC0*2*N*lPI 
60  TO  216 

219  0ISTT«YtI.J.K>-Yll.E8»NC8,N£8l 
216  1FINCS.C0.NI60  TO  220 

IF(N.GT.0IIP3(N<l)/2 

lFIN.lC.Q)lP*(N«2»/2 

0ISTZ»Z«I.J.KI-ZILCS»NC8,NCO>*.9*FL0AT« INCO-l»/2»-,9*FUOATCIP» 
60  TO  221 

220  OISTZmZ(I.J.K).ZILCO.MCU,NCai 

221  OIST2«01STX**2*01STY**2*OISTZ**2 
RETURN 

END 
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•001  i 
00021 
•0031 

•oooi 
ooos: 
•oo«i 
0007  i 
••••: 
oooo: 
00101 
•out 
0012: 
oois: 
0010: 
00131 
00141 

0017: 

ooio: 

001*1 

00201 

00211 

00221 

00231 

002*1 

00231 

00241 

00271 

00201 

002*1 

00301 

00311 

00321 

00331 

003*1 

00331 


SUBROUTINE  CONVERT U.M.N. LEO. NEO.NESI 
C0NN0N/SIIE/LX.LV,L*.NN.NY0T.L2«».L*SC0N 

CONNON/STRESS/OISTX.DISTT •OlBTZtDlSY  .OlSTt.Cl.RVSS.RVTY  «RVJ2 
C SUBROUTINE  CONVERTS  L.N.ANO  N INTO  EOUIVALENT  VALUES  FOR  RCRIOOIC 
C BOUNOART  CONDITIONS  ANO  FINOS  OlSTANCE  FROM  SCANNED  TO  TCSTEO  R ARTICLE 
C DETERMINE  LEO 

IF (L  .01.  1 .ANO.L.LC.  LSI  SO  TO  103 
IFU.LT.ll  00  TO  102 
LEB*M00<L«1«LXI*1 
00  TO  110 

102  INalMABSUI/LX 
LC0*l**LX»tA8S<L> 

00  TO  110 

103  LEB*L 

C DBTERN1NE  NEB 

110  IFlM.SC> l « AND.N.LC.LY  I BO  TO  HO 
IFIN.LT. II  00  TO  112 
NE0*N00<N»1.LYI«1 
SO  TO  120 

112  INalMABSINI/LV 
NEO*lN*LT«I ABS1NI 
SO  TO  120 
110  NEO*N 
C DETERMINE  NEB 

120  1F(N,sE.1>AN0.N.LE.L1«I  00  TO  120 
1FIN.lT.1I  00  TO  122 
NCB*NOO(N«l.L2*l*l*L2*»L2SCOM 
60  TO  130 

122  1N»1*IABS1NI/LI* 

NEB*1N*L<*«IABS1NI 
SO  TO  130 
120  NE8*N 
130  CONTINUE 
RETURN 
CNO 


38 


oooi : 

uoo2:  c 

0003:  c 

oooi: 

0009 : 

ooo« : 

ooot: 

oooo: 

oooi: 

ooio:  c 

oon: 

ooia: 

oois: 

oon: 

00191 

ooio:  c 

ooi?: 

ooio: 

ooio: 

0020: 

oooi: 

0022: 

0023: 

0021: 

0029: 

oo2« : 

0027: 

0020: 

0029: 

ooso: 

003i: 
0032: 
0033:  c 

oon : 

0039: 

0036: 

0037: 

ooso: 

0039: 

ooio: 

oon: 

0012: 

0013: 

oon: 

ooi9: 

ooio: 

0097: 

ooio: 
0019: 
0090: 
oooi:  c 

0092: 

0093: 

oo9i: 

0099: 


SUBROUTINE  0ISTFN|K1«K2'UV> 

calculates  DISTRIBUTION  FUNCTION  for  VX'  Vy,  and  vz  with  interval 
OV  IN  THE  RCGI0N  BETWEEN  ANO  INCLU0IN6  PL*NES  Kl  ANO  K 2. 

COMMOn/S IZE/LX  «L T .LZ.NP.NToT <LZ4  tLZSEGI 
CONNON/VEL/VX109.09. JOO).vT(Ol.04,jOO) .VZ(09i01,300> 
CONMON/OIST/NIRI 100)  tNlLUOO)  iNIMAX 
DO  1 jsl.100 
NILIU)«0 

1 NlRI J) *0 

calculate  total  number  OF  PARTICLES  between  Kl  AND  K2 

LKaK2-Kl«l 

NKIKlaLKiNP 

RNM1K2*FL0AT INK1K2 | 

WAITEI4.99) 

WRITE (6* 100IK1 <K2«NK1K2«0V 
•••  Vx  DISTRIBUTION  ••• 

WRITE (6. 101 1 
WRITC(6«t02) 

NINAXal 

LZ1*2*K1-1 

LZ2*2*K2 

00  2 k*LZ1 <L22 

DO  2 J*l.LT 

00  2 1 * 1 • LX 

V»VX 1 1 • JtK) 

2 CALL  SORTIV.OV) 

00  3 lal.NIMAX 

VR*>  9*OV*FLOAT (2*1-1 ) 

VLa,9*DV*FL0AT< 1-2*1) 

RNRaFLOAT (NIR I I ) ) /RNK1K2 
RNLsFLOAT(NILI 1 1 )/RNKlK2 

3 WRITE(6.103)VR.RNR.VL.RNL 
•••  VT  DISTRIBUTION  ••• 

00  1 1*1.100 

NIL  1 1 ) *0 
1 NIR(I)*Q 
NIMAXal 

00  9 k*L21 «LZ2 
DO  9 jal.LY 
DO  9 1*1. LX 
V*VY ( I . J.K) 

9 CALL  SORTIV.OV) 

WRITE!  6.101) 

WRITE|6il02) 

DO  6 lal.NIMAX 
VR*.9*0V*FL0AT 1 2*1-1 1 
VL*.9*OV*FLOAT 11-2*1) 

RNRaFLOAT (NIR 1 1 ) >/RNKlK2 
RNLBFLOATINILI I ) I/RNK1K2 

6 WRITE (6« 103 )VR«RNR.VL.RNL 
•*•  V2  DISTRIBUTION  •*• 

DO  7 1*1.100 

NIL! I )*0 

7 NIR 1 1 ) *0 
NIMAXal 
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0036:  DO  a K»L21.L22 

0037:  oo  a j*i.lt 

oosa:  oo  a i*i.lx 

0039:  V*V2 ( 1 . J.K  > 

0060:  a CALL  SORT ( V iOV ) 

0061:  WRITE t6« 1051 

0062:  WRITE  <6. 102 ) 

0063:  00  9 ISI.NIMAX 

0064:  VR*.3*DV*FLOAT«2M-l» 

0065:  VL*.3»0V4FL0AT ( 1*2*1 ( 

0066:  RNRaFLOATiNIKtn >/RNKlK2 

0067:  RNL»FL0ATINIUII>/RNK1K2 

006a:  9 WRITE16.103IVR.RNR.VL.RNL 

0069:  99  FORMAT  I /////2X . 30 AT A FROM  SUBROUTINE  OISTFnB* 

0070:  100  F0RMATI2X.80ISTRIBUTI0N  FUNCTION  FOR  PARTICLES  BEThEEN  PLANES8 . 13 

007i:  2.8  AND  8. I3/2X .8NUMBER  OF  pARTICLCS  IN  SAMPLE  s 8. 

0072:  3IS.3X.80V  * 8.FI3.6) 

0073:  101  FORMAT <2X. 8***  VX  DISTRIBUTION  •••8) 

0074:  102  FORMAT ( 3X.8VAV68. 3X . SN/NT0TS.4X.8VAVG8.3X.8N/NT0T8 I 

0073:  103  FORMAT I1X.2F7.4.3X.2F7.4) 

0076:  104  FORMAT C2X.8***  VT  DISTRIBUTION  ••*») 

0077:  103  FORMAT <2X.8*«*  V2  DISTRIBUTION  «**8) 

ooto:  return 

0079:  END 


00011  SUBROUT  INC  FORCC 

00021  C0NM0N/NISC/6RMMA .TAU.H.FX.FT.FZ.R.RC.CP. ICALC .CCHECK 

00031  C0MN0N/HISC2/F0RCEX.F0HCET.F0RCEZ 

00041  CONNON/STRCSS/OISTX.DISTT  <OlSTZ.DIST.Dt3T2.Cl.RVXX.PVYT.RVZZ 

00091  C2aClaa2 

00041  IFIC1. CO. 0.0)60  TO  1 

000T1  FORCC X a IC2*Cl)*OISTX/OIST 

00001  F0RCCTa<C2«Cl)*0ISTT/0IST 

0009!  F0RCCZa<C2>Cll*01STZ/0IST 

00101  60  TO  2 

0011!  1 FORCExaO.O 

0012!  FORCCtaO.O 

00131  FORCCZaO.O 

00141  2 FXaFX.FORCCX 

0019!  FYaF TaFORCCV 

0016!  FZaFZ.FORCCZ 

00171  RETURN 

0016!  CNO 


41 


0001 I SUBROUTINE  NABORS 

0002 1 C FOR  ANT  PARTICLE  (I.J.RI  SUBROUTINE  CALCULATES  I«L.J-M,  AND  K-N  WHERE 
0003:  CIL.M.N1  IS  A NEIGHBOR.  VALUES  ARE  STORED  IN  ARRATS  IX  ( I , J.  1 1 . IK  < I . J. 2 1 
00041  CITII.J.l)  etc.  WHERE  I HERE  DENOTES  UST.2ND.  3RD  ETC.  and  j the 

0003:  c particular  one  of  the  neighbors,  two  lxprcssions  are  used  to 
OOOt:  C CALCULATE  IX.it.  ANO  it  DEPENDING  ON  the  value  of  k of  the  scanned 
0007:  C PARTICLE.  ANO  THE  3RD  INDEX  DETERMINES  WHICH  EXPRESSION  IS  USED. 

0000:  COMMON/SC AHCH/NNAX.NMIN I NN.RCR I T.RCRIT2 

0009!  C0MM0N/NN80RS/JXIA,24,2> . Iyl6.2H.2t .I2IA.2W.2I ,NUM|A) 

OOIOI  DO  30  1*1. A 

00111  30  NUN  1 1 ) »0 

00121  DO  100  LL*1.21 

00131  DO  100  MM* 1 .21 

uom:  00  100  NNN*  1.21 

0013:  L*Ll-ll 

001 A • M*MM» 1 1 

0017:  N*NNN-11 

00  is:  I*L**2*M**2»N**2»L*N»L*M»M»N 

00191  IF  I I .GT.  NNAX)  GO  TO  100 

0020:  NUM II) aNUM ( 1 1 ♦ 1 

0021:  IRX2*L»M 

00221  IRT2*L*N 

0023:  IR22*M*N 

00241  R*NUM III 

0029:  IXlI.K.ll«IIRX2-l»IAes(M00lIRX2-I.2l n/2 

002a:  IXII.K.2(*(IRX2*I AOS (MODI 1RX2 • 2 1 1 1/2 

00271  IT(I.K.ll*(IRT2>l«IA8S(M00(IRr2>1.2>  M/2 

0020 : IT(I.K«2)*(IRT2*IAflS(  MOD  IIRT2.2MI/2 

00291  12 1 1 «K . 1 1»2*IR22»I ABS IMOOI I RT2 .2  M 

0030*.  12  ( I.K.2)*2*IR22>l«IA0S<M0n<  IRT2-1.2M 

0031!  100  CONTINUE 

0032:  return 


oooi : 
00021 
oooa: 
ooo« : 
ooos: 
000&: 
ooot: 
oooa: 
0009: 
oojo: 


SUBROUTINE  NRAN31 I Xl <X2 • I 1 

C0NH0N/NISC/6AMNAiTAU'H.FX.FT«F2'R'REtEPiICAkCtECHECK 
X3aSQRT(-2.Q*Al06(URAN31(t> I I 
X4a6t 2091 0330 72 *UR ANSI  1 1 ) 

X2aX3*SlN  ( X*t  I 
XlaX3*C03(X4» 

X2aX2/SORT (OANMAt 
XIaXl/SORT (6AMMAI 
RETURN 
ENO 


0001!  SUBROUTINE  OUTPUT 

0002!  C SUBROUTINE  PRINTS  SINOLC-PARTICLE  DISPLACEMENTS  ANO  VELOCITIES 

0003:  CONMON/SIZE/LX.LY .LZ»NP*NT0T«LZ4 .LZSEG4 

00041  CONNON/VEL/VXI09.09.300) • VY ( 04 .04. JOO ) . VZ ( 04.04 <300 ) 

OOOS:  COMMON/POSIT/XI04, 04. 3001 •V<04, 04. 3001. Z(04«04. 3001 

0006:  C0MM0N/MISC/6AMMA . TAU.H.FX  «FY  »FZ  «R»RE  «EP« ICALC .CCHECK 

0007:  HRITEIS.10I) 

OOOS:  SUNVXsO. 

0009:  SUNVTsO. 

0010:  SUMVZsO. 

OOli:  MRITE (6. 102  I 

0012:  00  9 KS1.LZ4 

0013:  00  9 J*l.LY 

0014:  00  9 1*1 .LX 

0019:  0X*X<I,J.K»-FL0AT(I»M.-.3*FL0AT<N0D»K.4>/2> 

00is:  OY*Y(I.J.KI-FLOAT( J>M.«.S*FL0AT<M00<K-1.2»  I 

0017:  OZ*Z<I.U.KI-0.9*FLOATIIK>l|/2) 

001s:  SUNVX»SUNVX*VXI I . J.K) 

0019:  SUMVYsSUMVY*VY ( I < J.K ) 

0020:  SUMVZ*SUMVZ>VZ  1 1 « J.K ) 

0021 : NPLANE*(K-1)/2*1 

0022:  9 WRITE<6<lQ3>NPLANE.I.J<K.OX<DY.OZ<VX(I.J.K).VV(I.J,Ki.VZ(I.J.K> 

0023:  NRITE»6.104>SUMVX,SUMVY.SUMVZ 

0024:  101  FORNAT ( /////2X . BOAT A FROM  SUBROUTINE  OUTPUTBI 

0029:  102  FORMAT I2X.8PLB.SX. BIS. 6X.Bj8.6X.SKB.9X.a0X8.12X. 80TB. 13X.B0ZB.13X 

0026:  2.8VX8.13X.8VT8.13X.BVZ81 

0027:  103  FORMAT (12.91 I4.3X|«6(E12.4.3X)) 

002S:  104  FORMAT  I IX .BSUMVX  • 8.E14.6.3X.8SUMVY  * B.C14.6.3X.8SUMVZ  ■ 8. £19. 

0029:  261 

0030:  RETURN 

0031 : ENO 
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ooois 

0002: 

0003: 

oooh: 

0003s 

0000: 

0007: 

0000: 

0009: 

0010: 

oons 

0012s 

00135 

ooi«: 
00135 
00165 
0017S 
00105 
0019S 
00205 
00215 
0022S 
00235 
00245 
00235 
00265 
00275 
00205 
00295 
0030S 
00315 
00325 
00335 
00345 
00335 
00365 
00375 
00305 
00395 
00405 
00415 
00425 
00435 
0044  5 
00435 
00465 
00475 
00405 
00495 
00305 
00315 
00325 


SUBROUTINE  PLANEAV 

C SUBROUTINE  CALCULATES  AVERAGE  VALUES  OF  POSITION  ANO  VELOCITY  FUR 
C PARTICLES  IN  A GIVEN  PLANE. 

CONMON/S12E/LX ,LV .L2.NP.NToT .LZ4.L2SE64 
COMMON/ VEL/VX I04.04.3005.VYI04.04.300) .VZl04.04.300) 
COMMON/POSIT/X(04.04.300) .7(04,04.300 ). 2 (04.04, 300) 

MR! TE(6«100) 

WRITE (6.101 ) 

SUMVXaO . 

SUMVVaO. 

SUMOXsO.O 

SUMDTaO.O 

SUMVZaO. 

00  30  KS1.LZ4.2 
OZAV«0. 

OYAVaO.O 

OXAVaQ.O 

VXAV«0. 

VTA V=0 • 

VZAVaO. 

00  20  Jal.LY 
DO  20  1*1 , LX 

02AVa02AV>2(I,J,K)*2(I.J.K4l)-,3aFL0AT( (K*l )/2)*,5*FL0AT(K/2) 
DXAVsoXAV«Xtl.J,K)+X(I, J.K*1)*2.0*FLOAT(I 1*2.0 
2*0 .3*FL0AT (MOO(K*l ,4 1/2 )*0,3*FL0AT(M00(K.4 )/2 ) 

DTAVaoTAV»T( I,J.K)*T(I.J.K«1 I*2.0*FLOAT( J)»2.0 
2-0.5*FLOAT(MOO(K.2))-0.3*FlOAT(MOO(K*1.2I ) 
VXAV*VXAV*VX)I»J»K)*VX(I»J,K*1) 

VYAVavYAV*VY(I«J.K)*VY(I«J,K*l) 

20  V2AV*V2AV*V2<I.V.K|*V2<Z.J.K*1> 

SUMVX*SUNVX*VXAV 
SUN V YaSUMV Y+ V Y A V 
SUMOXsSUMOX*DXAV 
SUMO VaSUMOY  40 V AV 
SUMVZaSUNVZ+VZAV 
02 AV«02AV/FLOAT ( NP ) 

VXAV*VXAV/FLOAT(NP) 

VYAV*VYAV/FLOAT (NP) 

V2 AVs V2 AV/FLOAT ( NP ) 

OXAVaOXAV/FLOAT (NP) 

OYAVaOY AV/FLOAT ( NP ) 

NPLANfa (R*l ) / 2*1 

30  WRITE(6.102)  NPLANE.OZAV.VZAV.VXAV.VYAV.DXAVtOTAV 
WR I TE ( 6 • 1 03 ) SUMV2 , SUMVX • SUM V T , SUMOX , SUMDT 
103  F0RMAT(//2X.aSUMSa.T24.E14.6,3x.El4.6.3X.El4,6.3x,E14.6.3X.E14.6) 

100  FORMAT (/////2X,aOATA  FROM  SUBROUTINE  PLANEAV8) 

101  FORMAT ( IX ,8PLa«9X .GOZAVa. 13X .OvZAVO .13X.aVXAVa.13X.OVYAVa.l3X, 
2B0XAV8 . 13X « BOV  A V8 ) 

102  FORMAT (1X,I3.3X,E14.6.3X.E14.6,3X.E14.6.3X,E14,6,3X,E14.6«SX« 
2E14.6) 

RETURN 

ENO 
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SUBROUTINE  POTFOHtNilSCANlI 
C0NN0N/*I3C2/F0RCER'F0RCET,F0RCE2 

C0NM0N/3TRCSS/0I3TX iD13TT *01312 tDI8T«0I9TI «C It PVRR«PVTT  *PV22 

C 0MN0N/NISC/8ANNA • TAU*H«FX(FT «F2« A«RE «EP « tCALC«ECHCCK 

OlSTaSOMTlOISTtl 

ClaEXPI *R* IRE*OISTal *0 1 1 

IF < ICALC  *NE« 1 • ANO.N.LT • 1 1 Cl«0.0 

IF(ISCAN1.C0.2>  80  TO  139 

CALL  FORCE 

|F(ISCAN1,NC.3>  80  TO  IA9 
PVRR»PVRR*OISTR*FORCER 
PVTT>PVTT«OISTT*FORCET 
PV22apv22«DlST2*FoRCE2 
139  C«C1 

IF  I ICALC  ,NE.1.AN0.N,lT.1»  cao.o 
CPatP.0.9«0'9*IC>i.0l**2 
1»»  CONTINUE 

return 

ENO 


SUBROUTINE  SAVE 

C SUBROUTINE  SAVES  CURRENT  VALUES  OF  POSITION  IN  VECTOR  RS  FOR  USE  IN 
C SUBROUTINE  SCAN 

C0NN0N/3IZE/LX«LV'LZ«NP(NT0T>L2lilZS£l0 
CONNON/POS1T/XIO«. M'SQO).  VISA, M.SOOl. 2(0*. 0*. SOU 
CONNON/POSITS/XS I M ' 0«  « S00  > « VS  I ••»  • 9«  < 300 ) • ZS  < •*  « 0*  , SOI  > 

DO  LOO  Ml.LZO 
00  100  J*1 *LV 
00  100  I»1.LX 
XS(l«J*K|sXII*J«K| 

TSU«J.Kt«va«J.Kt 
100  ZSII'JtK>>Z(I«J«K) 

return 

END 


0001 1 
ooooi 

00051 
0000 1 
00051 
00001 
0007 1 
0000! 
00091 
00101 
OOUt 

ooioi 

00151 
00101 
00151 
OOUt 
00171 
00101 
00191 
0000 1 
00011 
0000 1 
00051 
00001 
00051 
0000 1 
00071 
00001 
00091 
00501 
00511 
00501 
00551 
005OI 
00551 
00501 
00571 
00501 
00591 
00901 
00911 
0090 1 
00951 
00991 
00951 
00901 
00971 
00901 
00991 
00501 
00511 
OOStt 
00551 
00591 
00551 


SUBROUTINE  SCAN! I • J.K.ISCAnI. ISCANO) 

C 11  •***•"«•«*  THC  NEAREST  NCtOHBOftS  TO  PARTICLE 

r Ht*"E8T  RCUHBORS  A*t  ASSUMED  TO  It  TMOSC 

C RARTICLC^/j  »rS.'Kl«  **  ,8C*"E«E»  T*«C  rONcC  EXERTED  ON 

CONNON/STRC3S/DlSTX.OISTT»OlST|*OIST»OI8TO.Cl»R9XX.BlfvT.NM»X 

CONNON/RQSI TS/X 1 09 ,09 . 500 ) , 7 < 09 , 09 , *00 » .0 U9 . Ia!  I , ' 

COMMON/MI SC/OANNA.TAU.H.FX.FVi 72 « R,N[, CR • ICALCtCCHTCM 
CONNOX/UMCH/WMl  iNXINt  NN , ACM  j T . OCR  1 TO 
COMMON/NNBORS/tx 14.09.0) . It  14.09 .2) . 10 (A <09 <0 1 NuNi A i 

OIMCNSION  HU*MN|09,09.500I  ^ 

rx»o.o 

FT«0.0 

F0»0.0 

CR«0, 

RVXXae.O 
RVTVSQ.O 
7W*OiO 
RK*R00IR«9|Al 
00  TO  1S.4.7.0).KK 
5 1N5X»1 
1N5T«0 
IN50>0 
00  TO  9 
0 IN5X«1 
IN5T»1 
IN50«1 
00  TO  9 
7 lN5R*a 
1N5T«0 
1R50>0 
00  TO  9 
0 1N5X*0 
IN5T«» 

INSUl 
* CONTINUE 

IOlNNIN  .CO.  01  00  TO  151 
100  00  ISO  LU1.NNIN 
N1NAXaNUNM.1I 
00  199  NU1.N1NAX 
L>l«lX(Ll.Nl.IN5X) 

N*U*1T1L1 tNl • INST | 

NaK«U(L1.N1.IN521 

CALL  CONVCRTtL.N.N.LCO.NEO.NCO) 

CALL  CALOISTIL.N.N.LCO.NCO.NCo.I.J.ri 
tall  rotforin.iscanii 

199  CONTINUC 

150  CONTINUC 

151  CONTINUC 

iriNNAX.CO.NNIN)  00  TO  500 
iriI5CAN0.C0.0l  00  TO  900 
000  NUNNN(1.J.K)«0 
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00 


NN|NPlaNMIN«l 

00  390  Ll»NR»NPl.NNAX 

NiNAXaNUMUll 

00  300  M|at,MlNAX 

L»1*IX<L1.N1*IN3X> 

A«J*ITU1.NI*IN3Y» 

N«K«I*<L1.N1'IN32I 

CALL  C0N¥ERTtL.R.N«LE0.NE0,NEQ» 

CALL  CALOISTIL.M.N.LCO.MCO.NEQ. I < J.K) 

1FI0IST2  .6T.  RCRIT2  .OR.  0IST2  »LC.  i.OC-101  00  TO  909 
NUMNN < I . J . K» aNUMNN < I • J. K> ♦ t 
MRITC40)  L.R.N.LEO.NEO.NCO 
CALL  POTFORTN.ISCANIT 
300  CONTINUE 
390  CONTINUE 
00  TO  900 

000  NIU**NUMNN(I.J.K> 

IF |NI JK  .CO.  01  SO  TO  900 
00  090  L1«1.NIJK 
REAOIO)  L.N.N.LCQ.NCQ.NCQ 
CALL  CALOISTIL.N.N.LCO.NCO.NCO. I • J.KI 
CALL  POTFORIN.ISCANII 
090  CONTINUE 
900  CONTINUE 
RETURN 
END 
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oooit 

00021 

00031 

ooo*: 

0003: 

ooo*: 

ooo?: 

oooo: 

0009: 

ooio: 

oou: 

0012: 

0013: 

ooi*: 

oois: 

ooi*: 

00171 

ooio: 
0019: 
0020 : 
0021: 
00221 
0023: 
002*1 
00231 
002* : 
00271 
00201 
0029: 
00301 
00311 
00321 
00331 
003*1 
00331 
003*1 
00371 
00301 
00391 
00*01 
00*11 
00*21 
00*31 
00**1 
00*31 
00**1 
00*71 
00*01 
00*91 
00301 
0031 1 
00321 
00331 
003*1 
00331 


SUBROUTINE  SOLVE 

C SUBROUTINE  SOLVES  EQUATIONS  OF  NOTION  TO  OBTAIN  VALUES  OF  THE 

C POSITION  ANO  VELOCITY  AT  TINE  TAU  ♦ H OIVEN  THtIR  VALUES  AT  TINE 

C TAU.  THE  PROCEDURE  IS  VIA  A FOURTH. ORDER  RUNCE.KUTTA  NETHOO. 
CONNON/SIZE/LX.LY.LZ.NP.NTOT.LZR.LZSEB* 

COMMON/ VEL/VX ( 0* • 0* • 300  > . V Y ( 0* . 0* . 300 ) . VZ ( 0* .09 . 300 ) 
C0NM0N/P03IT/X 1 0*. 0*. 300) . V ( 0* .0* « 300 1 .21 0* . 0* .300 ) 

CONMON/POS I T S/ XS < 0* . 0* . 300 ) • TS 1 0* • 0* • 300 1 • ZS 1 0* 1 0* t >00 1 
C0NN0N/3E ARCH/NN AX . NNIN , NN , RCR I T . RCR I T2 
CONMON/NISC/OAKHA.TAU.H.FX ,FT  «FZ .R .RE «EP. 1CALC .ECHCCK 
CONNON/OUNNT/PX ( 0* , 0* , 300 ) ,PT 1 0* .0* , 300 1 .PZ I OR . 0* . 300 1 .PVX 1 0* . 0* • 
2300). PWI 0*. 0*. 300). PVZl 0*,0*. 300). Xt(0*.0*. 300). VI (0*. 04,300). 
3ZI 1 0* 1 0* .300) . VXI ( 0* .0* .300 1 • VTl <0*.0*.300 1 . VZ1 (0*.0*«300 ) 

LCVEL  2.PX.PT .PZ.PVX .PVT.PVZ.XI .T1 .ZI.VXI.VTI.VZI 
DIMENSION  FZKO*. 04.02) 

C SAVE  VALUES  OF  POSITION  ANO  VELOCITY  AT  BES1NNIN*  OP  INTERVAL  IN 

C VECTORS  RI  AND  VI  ANO  INITIALIZE  PBS. 

00  SO  1*1. LX 
00  30  J*1.LY 
00  30  N*1.LZ* 

XKI.J.KIaXtt.J.K) 

Ttll.J.K)aY(I.U.K) 

ZIII.J.N)aZ(I.J.K) 

VXI(I.J.K)mVXll.J.N) 

VYI(I.J.K)*VY(I.U,K) 

VZI ( I .J.K)*VZ ( I . J.K ) 

PX( I . J.K)aQ,0 
PY|I.J,K)*0.0 
PZI1.J.K)*0.0 
PVXI 1 .J.K) *0.0 
PVt(l.J.K)*0.0 

30  PV2I 1. J.K)*0.0 

C ENO  SAVE  ANO  INITIALIZATION  SEQUENCE. 

00  900  N*l»* 

C SAVE  CURRENT  VALUES  OF  POSITION  IN  NS. 

CALL  SAVE 

C BEBIN  SEQUENCE  FOR  S0LVIN6  EQUATIONS. 

REMINO  0 
00  200  Kal.LZ* 

DO  200  Jal.LY 
00  200  1*1. LX 
IP<N  ,NE.  1)  00  TO  *0 
CALL  SCANIX.J.K.1.1) 

00  TO  *1 

*0  CALL  SCANI I. J.K.1.2) 

C SCAN  RETURNS  VALUES  OF  FX.FY.  AND  FZ.  THE  COMPONENTS  OF  FORCES  ACTINB 

C ON  THE  PARTICLE  BE1N0  SCANNED. 

*1  FXa2.0*RC*R*FX 
FYa2.0*RC*R*FT 
FZa2.0«RE*R*FZ 

IF (K  .EB.  1 .ANO.  M .EQ.  1)  FZI 11 .0. 1 )*-FZ 
IF  IK  .EQ.  2 .ANO.  M .EQ.  1)  FZ! 1 1 • J,2)*«FZ 
IPIICALC  .NE.  1 .ANO.  K .EQ.  1)  FZaO.O 
IP ( l CALC  .NC.  1 .ANO.  K .EQ.  2)  FZaO.O 
CO  TO  (103.110. 113.120I.M 
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009* : 
0097: 
00901 
0099: 

oooo: 
oo*i : 
00*2: 
00*3: 
00*9: 
00*91 
oo**: 

00*7} 

0068: 

00*91 

0070: 

00711 

0072: 

0073: 

0079: 
00791 
007*: 
0077: 
0078: 
0079: 
0080: 
0081 : 
0082: 
0083: 
0089 1 

00*9: 

008*: 

00871 

00881 

00091 

oo9o: 

00911 

0092: 

00931 

00991 

0099: 

009*1 

00971 

00981 

00991 

0100: 

0101: 

0102: 

01031 

010*: 

01091 

010*: 

01071 

0100: 

0109: 

0110: 


109  PVXU, J<K1SPVX<I< J.KMFX/8,0 
PVT  11, J<KISPVT<I< J.KMFY/6.0 
PVZ<l.J.Kt8PVZ<I,J,KMFZ/*.0 
PX<I.J,KlsPX<I.J,KMVX<I,J,Kl/*,0 
pt<i,j,kispt<  i.j.kmvyu.j.k  1/6.0 

Pi<  I<J<KISPZ<  I.J.KM  VZ<  I.  J.K  I /*,0 
X<I, J<K1BXI<I< J.K MH8VX<I<J<K 1/2.0 
T<I<J<KISY  I <I<J<KMH8VY<  I,  J.K  1/2.0 
Z(I<J,K)s/I(I<J,K) *H*VZ II,J,K)/2,0 
VX1I<J<KI8VX1<1<J<KI *HaFX/2  <0 
VT <I<J<K1SVY l < I, J.K MH8FY/2.0 
VZ<I<J<KIXVZI11 < J.K l*H*FZ/2.0 
60  TO  200 

110  PVX<l.J,KlaPVX<I.J.KMFX/3.0 
PVY<I<J<KI8PVY< l . J.K MPT/!.© 
PVZ<1<J<KI aPVZ < I < J.K 1 «FZ/3, 0 
PX<I,J.KIsPX<  I,  J<KMVX<I<J<K  1/3,0 
PT<  1<J<K1SPY<  I.J.KM  VT  II,  J,K  1/3.0 
PZ<1<  J<K18PZ<1<  J.KMVZU.J.KI/S.O 
X<I< J,K 18XI<I<J<K  I ♦H«VX< I,J<KI/2<0 
T<  I, J<KI8TI<1<  J. KUH8VT<1<J<K  1/2.0 
Z< I. J,K1«ZI<I.J.K XH8VZ<I<J<K 1/2.0 
VX<l.J.K>8VXlll,J.K>9H*FX/2.0 
VT<1<J.K)xVY1(I.J.K)»H«FT/2.0 
VZ<I,J.KIsvZI<I,J,KMHsFZ/2.0 

60  TO  200 

119  PVX(I,J,KlsPVX<I,J,KMFX/3,0 
PVY<I,J.KlaPVY<I.J,KMFT/3.0 
PVZ<I,J.KlaPVZ<l.J.KMFZ/3.0 

PX< 1 , J, K |apx<!,J. KI9VX<1<J<K 1/3.0 
PT < I.  J.KlsPY <I<J<KMVT<  I,  J.K  1/3.0 
PZ<I <J,K 18PZ<1<J<K )»VZ(1,J,K 1/3.0 
X< I. J<KI8XI<I< J.K)«H* VX I l.J.K I 
Til, J.K)8T1<1.J,KI«H*VT<1,j,KI 
i < I, J, K I *ZI( I, J.K MH8VZ<1<J<K1 
VX<I,j.K|8VX1<1,J.K1«H*FX 
VT <I<J<K1SVVI <1<J<K18H8FY 
Vi < 1 . J,K I aVZI < 1 <J,KI*H*FZ 
60  TO  200 

120  PVX<I,J,K|#PVX<l,J,K!»FX/*,0 
PVT<I<J<KI aPVT  ll<J<Kl8FY/*.0 
PVi<I,J.KlsPVZ<I.J,KMFZ/*,0 

PX < 1. J, Klapxi I. J, K19VX<1<J<K )/*<0 
PT<  I.  J<K18PT<1<  J<KMVT<I<J<K  1/4,0 
PZ<  I.  J<K1SPZ<I<  J<KMVZ<I<J<K  1/6.0 
X<I<J<KI8X1<1<  J,KMH*PX<1,J.KI 
T<I<J<KI8TI<  1<J<KMH8PT<  I.  J.K I 
Z<  1,  J<K  I «ZI  <1<J<KMH8PZ<1<J<KI 
VX <1<J<K1SVXI<  I,  J<KMH8PVX<I<J<K I 
VT<  I,  J<KI8VT1<1<J<KMH8PVT  1 1,J,K1 
Vi  <1<J<KI8VZ1<  I <J<KMH8PVZ<1<J<K1 

200  CONTINUE 

900  CONTINUE 
REMINO  8 

C CALCULATE  WORK  OONE  IN  INTERVAL. 


CALL  SAVE 
MONK  aO.O 
00  *10  J«1.LT 
00  *10  1*1. L* 

CALL  SCANII.J. 1.1.21 
Fli«a,0«*C*MFl 

*19  N0RKaM0RK»0»9*tF2I<ItJ<l)»F2)*lltI«J«l>><I(ItJ«ll  I 
00  919  Jal tLT 
00  91*  1*1 (LA 
CALL  SCAN! 1<J<2<1<2I 
F2*-2.0*RE*R*F2 

91*  HOHK*MOAK«0.9*<F21<liJ.2MF2l*(2IItJ>2l>21II«Jt2n 
CCHCCK*CCHCCK+W0RK 

hcmino.a 

920  TAU*TAU«H 
RETURN 
ENO 
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0001 : SUBHOUT INC  SORT(V«OVI 

00021  C0NN0N/0IST/N1RI lOOttNlLllOO) «N1NAX 

00031  00  1 jsltlOO 

000*:  V1«0V*FL0AT I J»1 1 

0009:  V2*QV*FL0AT  t Jl 

000*:  IFIV1.LE.V.AN0.V.LT.V2I  60  TO  2 

ooot:  vin»dv»float t i* ji 

ooo*:  V2NbOV*FLOAT I • ji 

000*:  IF(V3N.Uf .V.ANO.V.lT.VIN)  *0  To  3 

ooio : t continue 

oou:  *0  to  « 

0*111  2 N1RIJ)bN1R< J>»1 

0013:  IFIJ.sT.NIHAXININAXaJ 

001*1  60  TO  • 

00131  3 NIt(J)«NIL(J)*l 

00l«:  ifij.bt.ninaxininaxbj 

00171  * CONTINUE 

001*:  RETUNN 

001*1  ENO 
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oooi: 

0002:  c 
0003:  c 
oooo : c 

00051 

ooo*: 

ooor: 

oooo: 

ooos:  c 

ooio: 

ooui 

0012: 

00131 

ooih: 

00131 

001*:  c 

ooir: 

ooio: 

0019: 

0020: 

0021: 

0022: 

0023: 

0020 : 

00231 
002*: 
00271 
00201 
00291 
00301 
0031 : 
00321 
0033! 
003*1  C 
0033! 
003*1 
00371 
00301 
00391 
00*0: 
00*11 
00*21 
00*3! 

00**: 

00*31 

00**: 

00*7: 

00*0: 

00*9: 

00301 

0031!  c 

00321 

00331 

003*1 

00331 


SUBROUTINE  STRUT  , tkifv til  POSITIONS  AND  ASSI6NS  THEN 

^rormo?fTitr*Cro;OINrTO°*,Ni;«tl:Ll*N  UISTRIBUTION  WITH  BERN  ZERO 
"no  S7RN0RR0  OEVUTION  ■ 1/SORT <S*NRR > ^ 

C0NN0N/SIZE/E*.LT.L2.NR.NT0T.U9*tZSE**  s00> 

CONNON/WEC/V*  < 0*  . 0*  . 300  . Vj < .z<ON.O« , 300 » 

CONNON/POSlT/X10H,0*.300».r«0*.0*«3OO  Z lc>cchcCk 

COMRON/NISC/GANNA«TAU.H.FX.FttFZ«H.Rl 

CALCULATE  INITIAL  POSITIONS 
00  7 k«1«L2* 

00  7 jal.LV 

^lIj^UFLORTin-l.O^O.S.FtORTlNOOlR.OT/a) 

Y«!.i:R»-FLORTU».1.0*0.3.FUO*T(ROO(K.1.2>» 

7 ZU*U.IO«0.3*FC°»T1<R»*M*» 

ASSIGN  INITIAL  VELOCITIES 

IRANsO 

kscrle>o 

SUMVXaO. 

SUNVTaO. 

SUNVZ*0« 

SUMV2aQ. 

DO  0 Ral.LZ* 

00  0 jal.LY 

00  * 1«1«LX  . . 

CALL  NRAN3UX1,X2, IRANI 

VX(I.J«R>»*1 

VT»l«J«KI,X2 

CALL  NRAN31 (XI «X2. IRANI 

VZlI.JtKI«Xl 

SUNVXaSUHVX*VX( l* J.KI 

SUNVT»SUNVT*VT1I*J«KI 

* w •«““  *>  "" 

0VX»»S«NVX/FL0AT(NT0T) 

0VT««SUNVT/FL0AT1NT0T1 
OVZ»*SUNVZ/FLOAT 1NT0T I 

C«l. 

9 CONTINUE 

00  10  K«1*LZS 
00  10  jal.Lt 
00  10  I*1»LX 

VX(IiJ«RI*C*VXII»J.H>  *DVX 
VY(I«J»RI*C*V7II • J.KI ♦OVT 
VZI 1. J.K>»C*VZ< I» J.Kl*DVZ 

Sl)NVXaSUNVX*VX  1 1 • J.R  I 
SUN*TaSUNVT*VT(I.J*RI 

„ 

' £i!ySK5m.*  "«« 

C**.*FL0AT(NT0TI/16ANNR*SUNV2I 

C«SORT<C) 

KSCALE*! 

OVX«0, 


UVJ»0. 

SUMVXaQ. 

SUNVVaO. 

SUMViaO. 

SUMVSaQ, 

00  TO  9 

II  ElNlTaSUMV3/a. 

Mftirctt'Ui 

NRITEI*«I3)EINIT  tSUMVX  tSUMvT  « SUN Vi 
13  F0RMAT<8X«B0ATA  FROM  SUBROUTINE  START*) 

IS  FORMAT (3Xi*INIT I AL.  VALUES  OF  KINETIC  CNERBY  AND  NET  VELOCITIES*! 
3AEI«,0) 

RETURN 

END 


OOOII 

00021 

0003! 

ooom 

ooos: 

00001 

coot: 

0000! 

0000! 

0010! 

00111 

0012! 

0013! 

0010! 

0013! 

0016! 

0017! 

0010! 

0019! 

0020! 

0021! 

0022! 

0023! 

0029! 

0023! 

002«! 

0027! 

0020! 

0029! 

0030! 

0031! 

00321 

0033! 

00391 

0033! 

OOSti 

0037! 

00301 

0039! 

0090! 

0091! 

0092! 

0093! 

0099! 

0093! 

00901 

0097! 

0090! 

0099! 

0030! 

00311 

0032! 

0033! 

0039! 

0033! 


SUBROUTINE  THERMO I K1 .K2.ET0T I 

C CALCULATES  AVERAGE  VALUES  of  Flo*  VARIABLES  BETMEEN  FLAMES  K1  ANO  K2 
COMMON/POSIT/* 109, 09.3001. 7(09,09. 300), 2! 09. 09. 300| 
COMMON/POSITS/XS(09.09.300|.T3(09,09.300>.2S«09. 09,300> 
COMMON/STRESS/OISTX .OISTT , OIST 2.0IST, 01 ST2, Cl .PVXX.PV7T.FV2Z 
COMMON/SIZE/LX.LY .L2.NR.NT0T .L29.L2SE69 
COMMON/ VEL/VX( 09, 09 .300! • vT 1 09 • 09 . 300 1 • V2 ( 09 ,09 , 300 1 
COMMON /MISC /GAMMA. TAU.H.FX ,F 7 .F2.R.RE.EP, tCALC.ECHECK 
C CALCULATION  OF  0ENSIT7  NORMALIZED  TO  2CR0*TENPERATURC  value. 

oens«o. 

2R1»0, 

2K2»0. 

L21*2*K1*1 
L22«2*K2 
00  1 jal.LV 
00  1 lal.LX 

2K182K1*2II.J,L21MZ(1,J,LZ1«1> 

1 2K2a2K2*Z ( I , J.LZ2 |«2( I , J«L22al t 
02a ( 2R2“2K1J /FLOAT (NP! 

Tla.SaFLOAT (K2-K1 ) 

0CNSaTl/02 

C CALCULATION  OF  AVERA6E  VELOCITV. 

VAVWO. 

NKlK2aNP* I K2*K1*1 I 
00  2 K8L21.L22 
00  2 jal.LT 
00  2 lal.LX 

2 VAVaVAV*V2 ! 1 , J.K ) 

VAVaVAV/FLOAT (NK1K2I 
ETRANSaFLOAT (NN1K2 |aVAVa*2/2. 

C CALCULATION  OF  TEMPERATURES. 

TXaO. 

TTao. 

TZaO. 

TaO, 

00  3 K8L21.L22 
00  3 J»1.LT 
00  3 lal.LX 
TXaTXaVX! I.J.K)**2 
T7aT7«VT(I,J.K)**2 
VZiaVZ ( I t J.KI-VAV 
• 3 T2aT2*V21««2 

E THERM" !TX*TT*T2!/2. 

TXaTX /FLOAT (NK1K2*3 ) 

TTbTT/FL0AT<NK1K2*3» 

T2=T2/FL0AT (NK1K2*SI 
TaTX»T7*T2 

C CALCULATE  KINETIC  CONTRIBUTION  TO  PRESSURE  TENSOR 
PKXXAVa3.0*DENS*TX 
PK7TAva3,0*OENS*TT 
PK22AVa3.0aOENS«T2 

C CALCULATE  AVERAGE  POTENTIAL  ENERGT  ANO  POTENTIAL  CONTRIBUTION  TO 
C PRESSURE  TENSOR 
EPAVao, 

PVXXAVaO. 


/ 


5o 


0036: 

OOSTI 

ooss: 

0059t 

oooo: 

oo6t: 

0062: 

00631 

0069: 

0063: 

0066: 

0067: 

006«: 

0069: 

ooro: 

007i: 

00721 

0073: 

00791 

0073: 

0076: 

0077: 

oo7a: 

00791 

oooo: 
ooai : 
0002: 
ooo3: 
0009: 
ooos: 
00061 
0007: 

oooo: 

0009: 

0090: 

00911 

0092: 

00931 

0099: 

0093: 


PVYYAVaU. 

PVZZAVaQ. 

CALL  SAVE 

00  9 k«L£1'L22 

00  9 jal.LY 

DO  9 |al,LX 

CALL  SCANII<J<K<3<11 

PVXXAVaPVXXAV*R«RC*0ENS*PVXX/FL0AT(NKlK2i 
PVTYAvaPVYYAV»R*RC*0CNS*PVYY/FL0AT)NKlK2l 
PV22AVaPVZZAV*R*RE*0ENSaPV2Z/FL0AT<NRlK2l 
9 CPAVatPAV»EP 
REMIND  • 

EPOTaLPAV 

ETOTaETOT ♦CTRANS»ETHERN<£PoT 

EPAVatPAV/FLOAT INK1K2I 

PXXAVaPKXXAV«PVXXAV 

PVYAV*PKYYAV<PVYYAV 

PZZAVaPKZZAV<PVZZAV 

MR1TEI6.99I 

WRITEI6<1001K1<K2<NK1K2 
WRITE) 6<101I0ENS<VAV<EPAV 
WRITE (6< 102 »ETRAN3,ETHERN,EP0T 
WRITE (6< 1031  TX<PKXXAV<PVXXAV<PXXAV 
WRITE )6< 109 ) TY<PKYYAV<PVYYRV<PYYAV 
WRITE (6< 105 ( T2<PK2ZAV<PV22AV<P22AV 
NRITE)6<1061  T 

99  FORMAT ) /////2X<S0ATA  FROM  SUBROUTINE  THERMO-VALUES  AVERA8E0  BCTWE 
2EN  PLANES  K1  AND  K281 

100  FORMAT) IH0<  8KI*8< 19  <2X<8K?aB< 19 < 3X  < 8NX1K2*0< 19  I 

101  FORMAT l3X<B0CNStTYaa<C13. 6. 3X<0V2AV0aW<C13.6<3X<0CPAV0a*<ElS,6» 

102  FORMAT (6X<BETRANSaa,E13.6<2X<BCTMERMaB<El3. 6<9X<BEP0TBB<E13<6I 

103  FORMAT ) 3X  <STXaSfE13.6<3X  <BPKXXaO<tl 3.6<3X<BPVXX«B<E13,6<  3X  tBPXXal 

2<E13.6I 

109  FORMAT 13X<BTYaS<E13.6<3X<BPKYTa8<El3.6<3X<8PVYYaB<C13.6<3X<WPYYaB 
2<E19«6> 

103  FORMAT  1 3X<BT2a8<E13.6<3x<8PK22a8<El3.6<3X<8PV22aa<C13<B<SX<BP2Za8 
2<El3<6> 

106  FORMAT  1 3X  <8TSB<E13<6I 
RETURN 
END 


57 


- ■ — - — - 


OOOli 

oooai 
0003: 
oooo : 
0003: 
oooo: 
ooot: 
oooo: 
0009: 
ooio: 
oou: 
ooxa: 
oois: 
ooio: 
oois: 


FUNCTION  UNAN31 1 1 1 
IFMUO'U.IO 
ii  laimuu 
10  J»I 
J*J*I3 

J»  J-  < J/07100000 ) *0710000)1 

J*J*J3 

J*J«(J/07100000 1*07100000 

J*J*S 

J* J- « J/07100000 ) *07100000 

Al»J 

I*J 

UMN3X*Al/07100000, 

NCTURN 

END 
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